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ABSTRACT 

We use statc-of-thc-art hydrodynamical simulations combined with a ID radiative 
transfer code to assess the extent to which the highly ionized regions observed close to 
z ~ 6 quasars, which we refer to as near-zones, can constrain the ionization state 
of the surrounding IGM. We find the appearance in Lya absorption of a quasar 
H II ionization front expanding into a neutral IGM can be very similar to a classical 
proximity zone, produced by the enhancement in ionizing flux close to a quasar em- 
bedded in a highly ionized IGM. The observed sizes of these highly ionized near-zones 
and their redshift evolution can be reproduced for a wide range of IGM neutral hy- 
drogen fractions for plausible values of the luminosity and lifetime of the quasars. The 
observed near-zone sizes at the highest observed redshifts are equally consistent with 
a significantly neutral and a highly ionized surrounding IGM. Stronger constraints on 
the IGM neutral hydrogen fraction can be obtained by considering the relative size 
of the near-zones in the Lya and Ly/3 regions of a quasar spectrum. A large sample 
of high quality quasar absorption spectra with accurate determinations of near-zone 
sizes and their redshift evolution in both the Lya and Ly/3 regions should confirm 
or exclude the possibility that the Universe is predominantly neutral at the highest 
observed redshifts. The width of the discrete absorption features in these near-zones 
will contain important additional information on the ionization state and the previous 
thermal history of the IGM at these redshifts. 

Key words: cosmology: theory - intergalactic medium - quasars: absorption lines - 
methods: numerical - radiative transfer - H ii regions. 



1 INTRODUCTION 

Unravelling the epoch of hydrogen reionization is one of the 
major remaining goals of observational cosmology. The dis- 
covery of extended high opacity regions in the absorption 
spectra of z > 6 quasars has provoked an intense discus- 
sion as to whether or not the tail end of the reionization 
epoch has been detected. Following the three year WMAP 
results, the complementary constraints on the start and 
duration of the reionization epoch from cosmic microwave 
background data are now actually very weak dPage et alJ 
I2OO6I : ISpergel et alj|200^ . An extended period of reioniza- 
tion starting at 2 ~ 15 and a rather rapid reionization 
around 2 ~ 6 both appear to be consistent with the data. 
Much hope for further progress is based on upcoming ex- 
periments to detect high-redshift 21 cm emission/absorption 
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JScott fc Reejll990l: iMadau et al1ll997l: IShaver et al]ll99gl : 
iTozzi et al.l l200d) . Currently, however, the stu dy of Ly- 
man series absorption in the spectra of quasars (•Fan et alJ 
2002, 2006b) an d garama-ray bursts fciardi & Loebl l2000l : 
Barkana fc Loebl 120041: iTotani et alJljoOQl are the key ob- 
servational probe of the hydrogen reionization epoch. The 
absence of a Gunn-Peterson trough jGunn fc PetersoJl965l) 
in quasar spectra at z < 5.5 shows unambiguously that the 
volume weighted neutral hydrogen fraction in the intergalac- 
tic medium (IGM) is less than one part in 10'* iFan et alJ 
I200I l2006b>l . Hydrogen reionization is therefore definitely 
complete at this time. However, at 2 ~ 6 the detection 
of dark absorption troughs in quasar spectra indicates the 
IGM neutral hy drog en fracti on is inc reasing rapidly with 
look back time (IDjo rgovski et al. 2001; Becker et aL[j20W|; 
Songaila fc Cowid 12002.: .Pentericci et al.ii2002l: 1 White et alJ 
200a:lFan et alJl2006bl) ' Whether a large fraction of the vol- 
ume is still neu tral at this time is, however, rather contro- 
versial. iFan et al.. („.2Q06bil obtain a lower limit on the vol- 
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ume weighted neutral hydrogen fraction of /hi > 10^'^'^ at 
z > G from a detailed analysis of the Lyman series absorp- 
tion in quasar spectra. There are also claims for an acceler- 
ation i n the upwar ds evolutio n of the Ly o; optical depth at 
z ~ 6 iCen fc McD onald 20o3: lFan erall [2002: Wh ite et alJ 
l2003l : I^i^^l.ll2006h) which may oc cur during the epoc h 
when individual H ii regions overlap (lGnedin|l2000l |2004|) . 
Other authors, however, argue that the data is consistent 
with a smoo th evolution due to the gradual thic kening of the 
Lya forest JSoneailal 12003: iBecker et ai]l2006l) . Large vari- 
ations in the transmitted flux along different lines-of-sight 
further complicate the picture. Large Lyman seri es troughs 
are de tected in several quasar spectra at 2; ~ 6 iFan et alJ 
|20063). Yet the spectrum of J1148 -|- 5251 at z = 6.42 ex- 
hibits patchy transmission, indic ating the IGM is highly 
ionized in at least some locations jWhite et alj|2003 . l2005t 
[Oh & Furlanetto 2005). This may be indicative of spatial 
fluctuations in the UV background which are expected at th e 
end of reionization dWvithe fc Loebl2006l:fFan et al.l2006bl) . 
but may also be attributed to t he fluctuations of the IGM 
density alone llLidz et al.1 120061: iLiu et alJl2006l) . It is thus 
currently unclear whether the IGM is highly ionized, sub- 
stantially neutral or indeed a combination of both just above 
z = 6, depending on the line-of-sight considered. Answering 
this question conclusively will have important implications 
for the nature of hydrogen reionization. 

The main limitation of Lyman series absorption as a 
probe of the IGM neutral hydrogen fraction is it ceases to 
be sensitive at /hi ^ 10~^ . Other approaches to probing the 
IGM neutral hydrogen fracti on at 2: > 6 include analysing 
metal absorption systems jOhl I2OO2I : iFurlanetto fc Loebl 
l2003l : iBecker et alj l2006lJ). the luminosity function of 
Lyg emitting galaxies jMalhotra fc RhoadjEooi : I Hu et alJ 
l2005l: lKa s hikawa et al."2006^. and the dark gap width distri- 
butioni^aschos fc Norman 2005; Gallcrani ct al. 2006) al- 
tho ugh all have their associated interpretative uncertainties 
(see iFan et alj|2006al for an recent review of results from all 
these techniques). 

In this work we shall concentrate on yet another probe 
of the neutral hydrogen fraction of the high-redshift IGM. 
Recently, claims that the IGM has a large neutral hydrogen 
fraction at 2; > 6 have been made based on the observed 
sizes of the small regions of transmitted flux observed in 
quasar abso rption spectra just blue-ward of the quasar's 
redshift LWvithe fc Loebl l2004l: iMesinger fc Hairnan .2004i: 
iMesineer et alJl2004l : IWvithe et alJl2005tlYu fc LJ2005^ . In 
the regions close to the quasar intergalactic hydrogen is ex- 
pected to be highly ionized. The typical sizes of these regions 
observed in Lya a bsorption at 2 ~ 6 are 5 — 10 proper Mpc 
jFan et al.ll2006bl) . Throughout this paper we shall refer to 
these regions as near-zones. This accounts for the possibility 
that these regions may be either the observational signature 
of quasar H 11 ionization fronts expanding into a significantly 
neutral IGM, or classical pro ximity zones due to enhanced 
ioni zing flux near t he qu asar jBaitlik et alJll988^ . 

IWvithe et all (l2005l) use the size of spectroscopically 
identified near-zones observed around seven quasars at 6 < 
z < 6.42 to infer /hi > 0.1. An alternative analysis by 
IMesineer fc HaimanI (|20^ uses the difference between the 
sizes of the near-zones identified in the Lya and Ly/3 re- 
gions of the spectrum of a single quasar at 2 = 6.28 to 
obtain /hi > 0.2. They attribute the slightly smaller size 



of the Lya near-zone to the prese nce of a damping wing 
from a significantly neutra l IGM iMiralda-Escude fc ReesI 
ll998l:lMiralda-Escudell99^ . If these results are correct, they 
would imply a very rapid increase in the neutral hydrogen 
fraction over a short redshift interval, from /hi ~ 10^* at 
z = 5.5 to /hi > 10~^ at 2 ~ 6.4. In contrast, by compar- 
ing the r elative sizes of Lya near-zones observed around 16 
quasars, iFan et alj J2006bl) conclude the IGM neutral hy- 
drogen fraction increases by a factor of ~ 14 assuming the 
near-zone sizes scale as (1 -I- z)~^f^^^^ over the same red- 
shift range. They find no strong evidence for a significantly 
neutral IGM at 2 = 6.4. 

These results appear to be in direct contradiction with 
each other. However, there are a number of uncertainties 
which may explain the origin of these different conclusions. 
Using the properties of spectroscopically identified near- 
zones to determine the IGM neutral hydrogen fraction re- 
quires assumptions about the ionizing luminosity and life- 
time of quasars, as well as the density distribution of the sur- 
rounding IGM. The clumpy nature of the IGM will produce 
large variations in near-zone sizes from one line-of-sight to 
the next, with the underdense regio ns in the IGM setting the 
observed sizes of these regions fe.o. lOh fc Furlan etto 2005]). 
Further uncertainties arise from the possibility of fainter ion- 
izing sources clustered around the quasar and the exact de- 
tails of the radiative transfer through the inhomogeneous 
IGM. The sizes of the spectroscopically identified near-zones 
themselves are also uncertain. They are determined from 
generally noisy, moderate or low resolution spectra and the 
systemic redshift of the quasar is usually not known to a high 
degree of accuracy. It is also not immediately clear whether 
the observed sizes of near-zones correspond to the full ex- 
tent of the regions impacted by the ionizing radiation of the 
quasars {i.e. the region behind the quasar H 11 ionization 
front), a key assumption when inferring the IGM neutral 
hydrogen fraction from the sizes of near-zones. 

In this paper we investigate these issues in some detail. 
We model highly ionized Lya and Ly/3 near-zones around 
high redshift quasars using a large hydrodynamical simula- 
tion combined with an accurate scheme for radiative transfer 
through an inhomogeneous IGM of primordial composition. 
We shall consider the information one may realistically in- 
fer about the ionization state of the IGM from the observed 
sizes of Lyman series near-zones. In section |5| we develop a 
simple analytical model for the sizes of the spectroscopically 
observed near-zones. A brief overview of our radiative trans- 
fer simulations is presented in section O and in section |1| we 
examine the dependence of Lya and Ly/3 near-zones sizes 
on the IGM neutral hydrogen fraction, quasar age and ion- 
izing photon production rate. W e consider the evol ution in 
Lya near-zone sizes observed bv lFan et al.l i2006lJ) in sec- 
tion |S] and present our conclusions in section [5] Details of 
our radiative transfer implementation and treatment of un- 
resolved self shielded clumps are found in the appendices. 

Throughout this paper we adopt the following set 
of cosmological parameters {Qm,Q,A,ilbh^ , h, as,n) = 
(0.26,0.74,0.024,0.72,0.85,0.95), consistent with the com- 
bined analysis o f three year WMA P and Lya forest data 
JViel et alJl2006l: ISeliak et al1l2006l). and a helium fraction 
by mass of F = 0.24 (e.g. lOlive fc Skillmanll2004ft . Unless 
otherwise stated all distances referred to in the text corre- 
spond to proper distances. 
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2 H II IONIZATION FRONTS AND THE 
PROXIMITY EFFECT 

2.1 The extent of the regions impacted by quasar 
ionizing radiation 

We firstly consider the case of an H ii ionization front (I- 
front) expanding into a significantly neutral IGM. Consider 
the idealised case of an observer positioned along the line of 
sight to a monochromatic ionizing source, isotropically emit- 
ting A'' Lyman limit photons per unit time into a homoge- 
neous, pure hydrogen IGM with a volume weighted neutral 
hydrogen fraction /hi = jihi/wh. If this observer could mea- 
sure the exact size of the H ii region around the source at 
various times t during the source lifetime, the H ii region 
expansion rate would be described by (see appendix A for 
details) 



dRio] 
dt 



N - 



i7ri?fo„QHiin| 



^TT-Rion/Hin-H 



(1) 



where i?ion is the distance of the H ii I-front from the 
source, auii is the recombination coefficient for ionized hy- 
drogen and nn is the proper hydrogen number density at 
a fixed baryonic density normalised by the cosmic mean, 
A = Ph/Pb, 



riH = 7.0 X 10"'' A 



(2) 



For small values of /hi the "I-front" obviously produces only 
a very small change in the ionized hydrogen fraction within 
a region of radius i?ion. 

Solving equation Q for a source with age tq yields. 



Rion — Rs 



1 



exp 



1/3 



where _Rs is the Stromgren radius, defined as 

1/3 



Rs = 



3N 



(3) 



(4) 



and tree ~ (whcvhii)"^ is the recombination timescale. If 
tq <^ /niircc the solution to equation Q may be rewritten 
as 



Rion 



4.2 



iV V''f tq 



(A/hi)^/^ 12 X 1057 s-i J V 107 yrs 



Mpc, 



(5) 



where we have adopted some fiducial values for a quasar 
embedded in a uniform IGM. 

Equation Q has been used extensively to model 
the si zes of H ii regions around high rcdshift quasars 
(e.g. IShap i ro fc G iroux 1987; Donahue fc ShuU 19 87j; 
Madau et alTTl999l ; j^adau & Rees 200g| ICen fc HairnarJ 
2000. : .White et al.l l200a iWvithe fc Loebl l2004l : lYu fcTuI 
20051: lYj2005l : IShapiro et alJl200d) . While equation Q is a 
reasonable approximation for the expansion rate of an H ii I- 
front embedded in a significantly neutral IGM, it does not 
predict the residual neutral hydrogen fraction behind the 
I-front. However, the sizes of Lya and Ly/3 near- zones ob- 
served in quasar spectra will depend on the value of the 
neutral hydrogen fraction behind the H ii I-front. The value 



of Rion can therefore be different to the size of the near-zone 
spectroscopically identified in Lyman series absorption. For 
small values of /hi, -Rion can be significantly larger than the 
observed near-zone size. 



2.2 A simple model for the observed sizes of 
near-zones in Lyman series absorption 

The Gunn-Peterson optical depth llGunn fc Petersonlll965^ 
resulting from the scattering of redshifted Lya photons by a 
uniform distribution of neutral hydrogen can be written as. 



r{z) 



(6) 



where = 4.48 x 10"^® cm^ is the scattering cross-section 
for Lya photons, H{z) is the Hubble parameter and c is the 
speed of light. 

Now consider the effect of a luminous quasar on the 
surrounding IGM behind its H ii I-front. The extent of the 
observed near-zone in the Lya region of the quasar absorp- 
tion spectrum, R^, will be set by the limiting Lya optical 
depth, Tjiin = ~ln_Fiini, at which the observer can clearly 
distinguish between the presence and absence of t ransmit- 
ted flux. However, as discussed bv lFan et ai]||2006bft . consis- 
tently defining the boundary of a near-zone is difficult, and is 
complicated by fluctuations in the IGM density distribution, 
the UV background and the low signal-to-noise and resolu- 
tion of the quasar spectrum. In addition, at lower redshifts 
the IGM has non-zero transmission outsi de of the near- zone, 
making the definition of Ra ambiguous. iFan et al.. 1.2006bl) 
define Ra as the extent of the region where Fu^n > 0.1 when 
the quasar spectrum has been smoothed to a resolution of 
20 A. The smoothing of the somewhat noisy spectra is very 
reasonable and a common definition is obviously essential 
for comparing observed spectra consistently. However the 
smoothing of the spectra can result in an underestimate of 
the estimated near-zone sizes if associated narrow transmis- 
sion peaks beyond the region of continuous transmission are 
obliterated, or in an overestimate if transmission from the 
IGM ionized solely by the UV background is included. We 
shall therefore adopt a definition of R^ more suited to a 
discussion of the general scaling laws for near-zone sizes in 
the absence of such practical difficulties, and will come back 
to a de tailed comparison using the definition of iFan et alJ 
j2006bD later. Here we do not apply any smoothing and de- 
fine the edge of the near-zone simply as corresponding to 
the last pixel at which the spectrum drops below the nor- 
malised flux detection limit of — 0.1, corresponding to 
an optical depth detection limit Tum = 2.3. Note that this 
definition will break down when transmission from the IGM 
ionized solely by the metagalactic UV background becomes 
common. 

Rearranging equation (|SJ and setting r — mm yields the 
neutral hydrogen fraction required to reproduce the optical 
depth detection limit, 

/iir ^ 5.4 X W'^Aiim (g) , (7) 

where Aum is the corresponding normalised baryon density 
as fraction of the cosmic mean and Za is the redshift of the 
Lya near-zone boundary. Assuming the hydrogen behind 
the quasar H ii I-front is highly ionized and in ionization 
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equilibrium, which is a good approximation in the vicinity 
of a luminous quasar, the ionization rate per hydrogen atom 
needed to reproduce the neutral hydrogen fraction at the 
detection limit is, 



rhn 
HT 



f lim 
J HI 



(8) 



where x = 1-158 if hydrogen and helium are fully ionized 
by the quasar and x = 1 if hydrogen only is ionized, and 
Q!Hii(ra) is the case A recombination coefficient evaluated 
at temperature Ta, which is the IGM temperature at Za- 
Converting this to a specific intensity gives 



JhS" ^ 1.34 X 10"^ 



2 X 104 K 



erg s 



Hz-\(9) 



assuming the quasar has power law spectrum with spectral 



index Qg = 1-5, % ~ 1-158 and ami 
evaluated at = 2 x lO"* K Ja; 



2-51 X 10~" cm^ s~ 



t al-ill99 7:l- The largest 
observable size of the Lya near-zone will then correspond to 
the distance R^^^ from the source beyond which the source 
ionizing radiation field drops below Jh™, 



T->max ^ 



J_ / Nhp 



rlim 
^HI 



1/2 



(10) 



where hp is Planck's constant and is the number of ioniz- 
ing photons emitted by the quasar per unit time. Combining 
equations Q and IIOI I finally yields. 



3.14 / N y^v 

Aiim V2 X 10^'' s-i y V 



2 X 104 K 



T lim 

2.3 



7 



-9/4 



Mpc.(ll) 



Once the near-zone size reaches R^^^, even if the H ii region 
around the quasar keeps growing, the residual neutral hy- 
drogen in the ionized region beyond _R™™ will be sufficient 
to halt the advance of the observed near-zone edge. Ob- 
viously the above definition for the maximum size of the 
near-zone obviously only makes sense if the amplitude of 
the metagalactic UV background above the Lyman limit is 
smaller than JjJ™. Also note the size of the observed near- 
zone in the Lya region will always be smaller than or equal 
to -Rion. The important point to take from this simple argu- 
ment is the different scaling of the two radii, Rion oc N^^^ 
and ii™™ (X iV^/^, with the ionizing luminosity, and that 
i?™™ is independent of the neutral hydrogen fraction of the 
IGM and the lifetime of the quasar. One may also derive a 
similar expression for the maximum size of the observable 
Ly/3 near-zone, 

iir.2.5(^) (^) R- (12) 

where T/j is the IGM temperature at the edge of the 
Ly/3 near-zone at 2/3. 

Equations IIH and 11211 should be considered as il- 
lustrative approximations only; radiative transfer through 
a realistic IGM density distribution is required to address 
this fully. A number of caveats should therefore be kept 
in mind- The upper limits for the near-zone sizes given by 



equations IIH and 1121 1 do not take into account the pos- 
sible efi'ect of a damping wing from a significantly neutral 
IGM, which reduces the size of the Lya near-zone relative to 
Ly/3 ( Mcsinger & Haiman 200^. The presence of the dense, 
self-shielded clumps responsible for observed Lyman-limit 
systems will further reduce the upper limit for the observed 
sizes of the near-zone compared to the estimates given by 
equations lllH and 11211 - Additionally, the overlying lower 
redshift Lya forest should reduce the observed size of the 
near-zone in the Ly/3 region. Note that for a predominantly 
neutral IGM, the sizes of the observed Lya and Ly/3 near- 
zones will both trace the position of the H 11 I-front and 
should both be given by -Rion, whereas if the IGM is highly 
ionized, the distinction between the near-zone and the gen- 
eral Lya forest will become highly ambiguous. 

However, assuming equation IIH holds, for Rion > 
RT'' {Rion > i?;?'"') the observed size of the Lya (Ly/3) 
near-zone becomes independent of the neutral fraction of 
the surrounding IGM. The same is true if Rion = -Rs- The 
timescale for the H 11 I-front to reach _R™™ in a homoge- 
neous IGM is approximately, 



t ~ 4.2 X 10'' A/hi 



Rn 



3.14 Mpc 



N 



2 X 105'' s-i 



yrs 



(13) 



For a realistic clumpy IGM in the vicinity of a quasar this 
timescale should be longer by a factor of a few, but for 
/hi < 0.1 this is still considera bly shorter than the canoni- 
cal quasar lifetime of ~ 10^ yrs llMartini2004l : lHopkins et alJ 
E005.') . The time scale for the H 11 I-front to reach 7?™'"' is 
~ 16 times longer again. The size of the Ly/3 near-zone is 
therefore potentially much more sensitive to the IGM neu- 
tral hydrogen fraction- A similar argument applies for the 
scaling of the observed near-zone sizes with quasar age. 



2.3 Comparison of the analytical model to 
observational data 

The solid curves in both panels of figure show R^^^ and 
its dependence on the number of ionizing photons emitted 
per second by the quasar. A*', using our fiducial assumptions 
of [riim, as, T^, Za] = [2.3, 1.5, 2 x lO" K, 6] in equation llTTll . 
The three curves shown are for Aum = [0.4,0.7,1]. In the 
left panel the series of dotted and dashed curves show Rion 
computed using equation @ for /hi = [1,0.5,0.25,0.1] and 
a quasar lifetime of tq = 10^ yrs. In the right panel -Rion is 
shown for — [0.1, 0-25, 0.5, 1] x 10^ yrs and a fixed neutral 
fraction of /hi = 0.1. When computing -Rion, A = 1 has been 
assumed, although in reality A should be slightly larger on 
average due to the clumpy nature of the IGM. However, 
we find this does not make any significant difference to our 
conclusions, as we shall discuss later. 

The filled grey circles in figure represent the observed 
sizes of Lya near -zones around qu asars at 5.81 < z < 6.42 
as determined bv lFan et alj |2006b). In the instances where 
no Mg II redshift was available, £an ct al. (2006b) have ap- 
plied a correction to the quasar redshifts of Az = -1-0.02 
in order to take into account the systematic redshift off- 
set of high ioniza t ion line s. Note further that, as discussed 
above, iFan et al.l (l2006d) have smoothed the Lya spectra 
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Figure 1. Comparison of the observed sizes of Lya near-zones as determined by Fan et al. {2006b) (filled grey circles) and Wyithe et 
al. (2005) (open circles) to the two analytical models discussed in section \^ which predict the extent of the H II I-front and the size of 
the Lya near-zone at 2 = 6, respectively. The sizes are plotted as a function of the number of ionizing photons emitted by the quasar 
per second, N. The solid curves in both panels give R^^^ for Aiii^ = [0.4,0.7, 1.0]. Left: The series of dotted and dashed curves give -Rion 
assuming a uniform IGM with A = 1, a quasar age of tq = 10^ yrs and an IGM neutral hydrogen fraction of /hi = [1, 0.5, 0.25, 0.1]. Right: 
The dotted and dashed curves give -Rion assuming a uniform IGM with A = 1, an IGM neutral hydrogen fraction of /hi = 0.1 and quasar 
ages of tQ = [0.1, 0.25, 0.5, 1] x lO'^ yrs. 



to a resolution of 20 A before determining the near-zone 
sizes. The open circles with error bars are Lya near-zone 
sizes as determined bv lWvithe et alJ (|200^ around a sub- 
sample of the same quasars at 6 < z < 6.42, adapted from 
their table 1. IWvithe et al.l i2005l) " used the Lya near-zone 
sizes as a proxy for the extent of the quasar H ii I-front, 
and determine the extent of the near-zones to be at the 
redshifts where the transmission extending blueward of the 
Lya emission line becomes comparable to the noise. The 
size of Lya near-zone around J 1030 -I- 0524 they quote was 
increased to correspond to the size of the slightly larger 
Ly/3 near-zone observed around this quasar. We have in- 
stead adopted Ra = 4.55 ± 1.65, with an error bar based on 
an uncertainty of A z = 0.03 in the redshift of J1030 -I- 0524 
JWhite et al.l2003l). We h ave also omitted J1048-I-4637 from 
the ^Wv^^^e^L^^2005^ sample, as this is a BAL quasar 
which makes a meaningful measurement of the Lya near- 
zone size difficult (X.Fan, private communication). For all 
except one of the quasars (J1630 -|- 4 012, the lowest lumi- 
nosity quasar) the lWvithe et all ^2005^ v alues are consistent 
with the corresponding determinations bv lFan et all J2006bl) 
within the error bars shown. As discussed bv iFaj^^D 
^2006b^ and later within this paper, the observed near-zone 
sizes evolve rather ra pidly with redshift. The larger near- 
zone sizes within the iFan et al 1 J2006bl) sample which are 
omitted in the IWvith^^iLn2005F data are for quasars at 
z < Q. For both data sets, the conversion from ultraviolet 
AB absolute magnitudes, M1450, to A'^ has been made using 
a generic broken power law spectral energy distribution. 



-0.5 



u-" " (1050 < A < 1450 A), 
(A < 1050 A). 



(14) 



All photons above 1 Ryd are considered to be ionizing pho- 
tons. We have not plotted error bars for A'^, although the 
values will be somewhat unc ertain dependin g on the exact 
spectrum of the quasars (e.g. IYu fc LullioOSl) . 

As already discussed, the observed size of the near-zone 
in the Lya region should be the smaller of ii™'"' or i?ion. 
Clearly from figure our simple analytical argument sug- 
gests that the sizes of the observed Lya near-zones are con- 
sistent with both the expected extent of classical proxim- 
ity zones and an H 11 I-front expanding into an IGM with 
a significant neutral hydrogen fraction, for reasonable as- 
sumptions for the quasar ionizing luminosity and lifetime. 
The fact that equation ^ and llll l give similar values 
for iZion and R^'^^ is a coincidence which explains the dis- 
crepant interpretation of the observed sizes of near-zones 
in the literature. Taking a combination of 0.1 < /hi < 1 
and 10" < tQ < 10^ yrs results in an H 11 region simi- 
lar in size to the observed Lya near-zones, whereas equa- 
tion mil reproduces the observed sizes even for small /hi or 
large tq. For example, one can reasonably adopt values of 
/hi = 10""^ and t = 5 x lO' yrs, wit hin the observationally 
determined limits o f /m > lO"^ '^ jFan et al. 2006b) and 
10® < iq < 10* yrs (lMartinill20oi) . For a quasar luminosity 
of iV = 2 X 10^^ s~^, assuming Aum = 0.4 and Za = 6, ac- 
cording to equation llll l one expects to see a Lya near-zone 
of Ra — 7.9 Mpc, consistent with observed sizes. However, 
equation © predicts a value of 7?ion — 22 Mpc, much larger 
than any observed Lya near-zone, although in practise this 
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will be slightly smaller due to the clumpy nature of the IGM. 
If the size of the observed near-zone is nevertheless identi- 
fied with _Rion, the inferred IGM neutral hydrogen fraction 
is erroneously large. 

Therefore, the size of individual near-zones in the 
Lya region alone appears to offer little in terms of constrain- 
ing the ionization state of the hydrogen in the surrounding 
IGM. This will be confirmed by our detailed numerical sim- 
ulations including radiative transfer in section 2] However, 
as we will show, this situation improves if both the Lya and 
Ly/3 near-zones are considered in tandem. 



Table 1. Hydrodynamical simulation resolution properties. 
Listed in each column are the total gas particle number, the mass 
resolution per gas particle and the gravitational softening length. 



Gas particle 


Gas particle 


Softening length 


number 


mass [M0//1] 


[comoving kpc/h] 


120^ 


1.61 X 10^ 


10 


240^ 


2.01 X 10* 


5 


4003 


4.34 X 10^ 


3 



3 CREATING SYNTHETIC ABSORPTION 
SPECTRA INCLUDING NEAR-ZONES 

3.1 Radiative transfer implementation 

Our ID radiative transfer implementation is based on an 
updated version of the multi-frequency photon conserving 
algorithm used by iBolton et alJ (l2004l) . The advantage of 
this single ray-tracing approach is that we may run many 
radiative transfer simulations with only modest computa- 
tional resources. In addition, we achieve the high spatial and 
temporal resolution which is required t o accurately deter- 
mine post-reio nization gas temperatures iBolton et al.l20oS 
iTittlev fc Me iksin 2006). Technical details and the relevant 
resolution tests can be found in appendices B and C, to 
which we refer the interested reader. 



3.2 Hydrodynamical simulations 

We combine our radiative transfer implementation with 
cosmological density distributions drawn from a large 
hydrodynamical simulation. The hydrodynamical simula- 
tion was run u sing the parallel TreeSPH code GADGET-2 
JSpringeJiioOsI) . The simulation volume is a periodic box 
60/i~^ comoving Mpc in length and contains 2 x 400^ 
gas and dark matter particles. Star formatio n is included 
using the multi-phase model of Spring cl fc Hernauisli 
i2003r with winds disabled. The simulation was started 
at z = 99, with i nitial conditions g enerat ed using the 
transfer function of lEisenstein fc Hiil il999() and cosmo- 
logical parameters consistent with the combine d analysis 
of th r ee year WM A P an d Lya forest data iViel et alJ 
I2OO6I: ISeliak et all l2006l) . (f7m, !^a, !^2b/^^ /i, cts, n) = 
(0.26,0.74,0.024,0.72,0.85,0.95). The initial gas tempera- 
ture was set to T = 272.8 K and the gravitational softening 
length is 3/i~^ comoving kpc. To check the numerical 
convergence of our results we also run two further hydro- 
dynamical simulations with lower mass resolution. The 
different resolution parameters for all three simulations are 
listed in table all other aspects of the lower resolution 
simulations are identical to the 400^ run. 

A friends-of-friends halo finding algorithm was used to 
identify the 5 most massive haloes within the simulation vol- 
ume in three different snapshots ai z — [6.25, 6.0, 5.75]. The 
total masses (dark matter, gas and stars) and coordinates 
of the haloes in the z = 6 snapshot are listed in table |5| 
These halo masses may be slightly smaller than, or similar 
to, the expected mass of bright quasar host haloes at z ~ 6. 
The typical comoving space density of quasars at 2 ~ 6 



Table 2. Coordinates and masses of the five most massive haloes 
in the 400'^ simulation volume at ^ = 6. The halo coordinates are 
given in units of comoving kpc//i. 



^pos 


Vpos 




Halo mass [Mq /h] 


51030.5 


46853.1 


18807.6 


2.65 X 10^2 


55886.6 


49077.1 


12091.1 


1.42 X 10^2 


52392.1 


48011.6 


21267.5 


1.33 X I012 


52061.7 


47983.3 


22023.3 


1.25 X 10^2 


49415.3 


47008.2 


18859.4 


1.05 X 10^2 



is ~ 10"^ Mpc"3 jFan et alJl200ll) . which coincides with 
the expecte d abundance of dar k matter haloes with mass 
~ 1O"M0 iMo fo WhitFll2nn2h . However, recent observa- 
tions of molecular gas surrounding J1148-I-525 1 indicate its 
host halo mass may be closer to ~ IO'^^Mq iWalter et alJ 
[2904). Theoretical arguments also suggest that the quasar 
occupation fraction of massive haloes may be significantly 
less than unit y, perhaps favouring hos t haloes with lower av- 
erage masses ijVolonteri fc R,ee3l2006ll . In any case, we have 
checked that our results depend only weakly on the halo 
mass by running simulations centred on haloes of smaller 
mass. 

At each redshift, five lines-of-sight where extracted 
around these haloes, giving 15 different lines-of-sight in to- 
tal. To construct the density distributions for the radiative 
transfer runs we splice together separate lines-of-sight drawn 
randomly from the simulation volume. The first 30h~^ co- 
moving Mpc of a spliced line-of-sight is constructed using 
the density distribution drawn around one of the massive 
haloes, in which we assume the quasar resides. The remain- 
ing 120/i~^ comoving Mpc of the sight-line is made from two 
further sight-lines drawn randomly from the simulation box 
at the same redshift. This gives a total length of 150/i~^ 
comoving Mpc, which is large enough to easily accommo- 
date the observed sizes of Lya and Ly/3 near-zones at 2 ~ 6. 
Each line-of-sight consists of 3000 pixels with a size of 9.9 
kpc at 2; = 6. This provides a good degree of convergence 
for our radiative transfer simulations (see appendix C for 
details). For reference, the redshift extent of each line-of- 
sight at 2 = [6.25, 6, 5.75] is Az = [0.50, 0.47, 0.45]. Since the 
density distributions are drawn from a snapshot at a single 
redshift, the gas density is rescaled by a factor of (1 -f 2)"^ 
along each spliced line-of-sight. 
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Table 3. Spatially uniform UV background models used for 
setting the neutral hydrogen and helium fractions of the sur- 
rounding IGM. The first two columns list the specific in- 
tensity and ionization rate at the Lyman limit in units of 
10~^^ erg s~^ cm~^ sr~^ Hz~^ and 10~^^ s~^ respectively. 
The last column lists the resulting neutral hydrogen fraction for 
mean cosmic density at 2 = 6 assuming ionization equilibrium at 
T = 2 X 10^ K. A power law spectrum with a\y = 1.5 has been 
adopted. 



Model Jbj^ r^ia /hi 


















1 




1 


2.0 


X 10" 


-5 


5.6 X 10" 


-5 


3.6 X IQ- 


-1 


2 


1.0 


X 10~ 


.4 


2.8 X 10" 


-4 


7.3 X 10- 


-2 


3 


3.0 


X 10" 


.4 


8.5 X 10" 


-4 


2.4 X IQ- 


-2 


4 


1.0 


X 10" 


-3 


2.8 X 10" 


■3 


7.3 X 10- 


-3 


5 


3.0 


X 10~ 


-3 


8.5 X 10- 


'3 


2.4 X 10- 


-3 


6 


1.0 


X 10" 


-2 


2.8 X 10" 


-2 


7.3 X 10- 


-4 


7 


3.0 


X 10" 


-2 


8.5 X 10- 


-2 


2.4 X 10- 


-4 


8 


6.0 


X 10" 


-2 


1.7 X 10- 


-1 


1.2 X IQ- 


-5 


9 


1.0 


X 10- 


-1 


2.8 X 10- 


-1 


7.3 X 10- 


-5 


10 


2.0 


X 10" 


-1 


5.6 X 10" 


-1 


3.6 X 10- 


-5 



3.3 Initial conditions 

We compute the transfer of radiation through an inhomoge- 
neous hydrogen and helium IGM from a source which emits 
A'^ photons per second above the H i ionization threshold 
frequency, where 

•I I^HI 

The source is assumed to have a power law spectrum with 
index Os, typical of quasars, normalised by Lhi at the 
H I ionization threshold frequency, vsi- We adopt Os = 1.5 
for this study. The source is placed 10/i~^ comoving kpc from 
the edge of the line-of-sight; all hydrogen and helium within 
this radius is assumed to be fully ionized. The distribution of 
normalised densities A = ph/ph is kept constant during the 
radiative transfer simulations, while the gas density evolves 
as (1 -|- z)'^ as the ionization front progresses along the line- 
of-sight. This should be a reasonable approximation over a 
typical quasar lifetime. 

We vary the lifetime of the source, Iq, the number of 
ionizing photons emitted per second by the source, A^, and 
the neutral hydrogen (and helium) fraction of the IGM. To 
set the neutral fractions in the IGM for each run, we assume 
the gas is in ionization equilibrium with a spatially uniform 
UV background. The different UV background models we 
use are listed in table |H] along with the resulting neutral 
hydrogen fraction at mean density at z = 6. Note, however, 
that this does not necessarily correspond to the resulting 
volume weighted neutral fraction in a simulated sight-line, 
as on average these will be moderately overdense and have 
varying gas temperatures, although it gives a rough guide 
to the expected neutral hydrogen fraction for a given UV 
background. All UV background models assume a spectral 
index of Qb = 1.5, although the background spectrum may 
well be softer if it is dominated by galaxies. This assumption 
does not have a significant effect on the results of this paper 
apart from small changes in the neutral helium fraction and 



the temperature of the IGM. There should also be signif- 
icant spatial fluctuations in the UV ba ckground at z ~ 6 
dWvithe fc Loebll200el : lFan et al.l l2006b') but these will only 
be important at the edge of the quasar near-zone itself and 
only if the flux from the quasar and the metagalactic ioniz- 
ing background are similar there. In the latter case spatial 
fluctuations of the metagalactic ionizing background may 
affect the measured near-zone sizes. 

3.4 Constructing synthetic absorption spectra 

We construct synthetic spectra from the output of the ra- 
diative transfer simulations as follows. Each line-of-sight is 
rebinned to have A'' = 4096 pixels of proper width 5R, each 
of which has a neutral hydrogen number density riHi, tem- 
perature T, peculiar velocity Upec and Hubble velocity 
associated with it. The Lya optical depth in each pixel is 
computed assuming a Voigt line profile, such that the opti- 
cal depth in pixel i corresponding to Hubble velocity WH(i) 
is 

JV 

Here h = {2k-BT /m-n)^^^ is the Doppler parameter, ctq = 
4.48 X 10~^* cm^ is the Lya s cattering cross-section and 
H(a, x) is the Hjerting function iHiertinelll93Sl 

a f°° e-y^ 
H{a,x) = - <^y^ (16) 

where x = [v-nii) - u{j)]/b{j), u{j) = whO') -I- Vpadj), 
a = AaAa/47r6(j), Ac = 6.265 x 10* s~^ is the damping con- 
stant and Aq, = 1215.67 A is the wavelength of the Lya tran- 
sition. For Ly/3 optical depths we adopt the same procedure, 
using instead = 7.18 x 10"^^ cm^, A/j = 1.897 x 10** s"^ 
and Xjs = 1025.72 A. In addition, for the Ly/3 spectra we 
also add the optical depth contribution from the Lya for- 
est at lower redshift, z' — [A;3(l -I- z)/\a] — 1, where z is 
the redshift of the Ly/3 pixels. The total Ly/3 optical depth 
is then r^°'(z) — tis{z) -\- Ta{z'). In order to model the 
Lya forest at z' , we use the standard technique of rescal- 
ing the optical depth distribution of synthetic spectra, also 
drawn from our 400"^ simulation at lower redshift, to re- 
produce the observed mean fiux (F) of the Lya forest (see 
Bolton et al. 2005 for details). The values f or the mean fiu x 
we adopt in this paper are interpolated from lSongailal i2004l) : 
(F) = [0.149,0.210,0.223] at z' = [5.12,4.91,4.70], corre- 
sponding to z = [6.25, 6.0, 5.75]. 

For comparison to observational data, we must also 
mimic the typical resolution and noise of quasar spectra at 
z ~ 6. The values for these vary in the literature from one 
line-of-sight to the next. We adopt a spectral resolution of 
4 A per pixel by rebinning our simulated spectra and adding 
Gaussian distributed noise with a total sig nal-to-noise ratio 
of 20 per pixel at the continuum level (e.g. iFan et al.l.2000'1 
and a constant read-out signal-to-noise of 80 per pixel. Ex- 
ample synthetic Lya and Ly/3 spectra are shown in figure |5| 

3.5 Hydrodynamical simulation resolution 

To simulate the propagation of ionizing radiation through an 
inhomogeneous IGM, the structure of the IGM must be suf- 
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Figure 2. Examples of the synthetic Lya and Ly/3 near-zone spectra used in our study. The spectra are drawn from the line-of-sight data 
shown in figure|^ The quasar is situated on the right of the diagram, and has an age = lO'^ yrs and ionizing luminosity N = 2x lO^'^ s~^ . 
The surrounding IGM neutral hydrogen fraction is set using UV background model 4, which gives a volume weighted value of /hi ~ 10~^. 
Here the Lya spectrum has been rebinned to a resolution of ~ 4 A per pixel and Gaussian distributed noise with a total signal to noise 
of 20 per pixel has been added, typical of moderate resolution quasar spectra. The Ly/3 spectrum has been rebinned to have the same bin 
size as the Lyo spectrum in velocity space. The solid line in each panel marks the extent of the quasar H II region, some way ahead of the 
edge of spectroscopically identified Lya and Ly/3 near-zones. These both lie at a distance of ~ 8.3 proper Mpc from the quasar. The dotted 
line shows the adopted flux threshold, = 0.1, used to identify the edge of the near-zone. 




54321 54321 54321 

[proper Mpc] [proper Mpc] [proper Mpc] 



Figure 3. Left: Density distributions at z = 6 drawn from hydrodynamical simulations with differing mass resolution. The dotted, dashed 
and solid lines correspond the 120'^, 240'^ and 400^ simulations respectively. The corresponding mass resolutions are listed in tabled Centre: 
The neutral hydrogen number density normalised by the total hydrogen number density as predicted by the radiative transfer algorithm 
for a quasar within the massive halo situated on the right hand side of the density distribution. Right: The resulting Lya near-zone profile. 
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ficiently resolved to correctly model the numb er of ionizing 
photon absorption event s per unit volume (e. q . iKohler et alJ 
l200fll: llliev et al.ll2006bl: lOnedin fc Fanll200fil^ . We test this 
using the hydrodynamical simulations listed in tabled We 
firstly select one of our spliced lines-of-sight at 2 = 6 from 
the 400^ simulation, and construct lines-of-sight from the 
same location in the lower resolution simulations. For all 
three lines-of-sight we then use the radiative transfer code to 
model the propagation of ionizing radiation from a quasar 
with iV = 5 X 10^'' s~^ and Iq = lO' yrs into fully neu- 
tral hydrogen and helium gas. We subsequently extract the 
Lya absorption spectra. 

The results of this test are shown in figure |3] From left 
to right the panels show the density distribution, neutral 
hydrogen fraction and the resulting Lya near-zone profile 
drawn from the radiative transfer runs in all three sight- 
lines. The degree of convergence in the H ii I-front posi- 
tion is good. The highest resolution simulation resolves a 
greater number of overdense peaks, but the effect that this 
has on the I-front position is counteracted by better resolved 
voids. The main difference is that one sees more detailed 
structure in the higher resolution simulation spectra. Note 
also the smooth nature of the Lya near-zone profile due to 
the damping wing from the neutral IGM; the edge of the 
Lyg near- zone lies behind the H ii I-front as a result {e.g. 
iMesineer fc HaimanI [20041) . We conclude that the 400^ hy- 
drodynamical simulation is sufficiently resolved for the pur- 
pose of this work with one caveat; the very dense regions 
responsible for Lyman limit absorption systems will not be 
sufficiently resolved. We will come back to this in section 5. 



4 MULTI-FREQUENCY RADIATIVE 

TRANSFER THROUGH COSMOLOGICAL 
DENSITY FIELDS 

4.1 An example line-of-sight 

In figure |21 we had shown synthetic absorption spectra for 
an example line-of-sight at z = 6. Figure |1| shows a variety 
of physical properties along the same line-of-sight together 
with a high resolution, high signal-to-noise version of the 
Lya region of the absorption spectrum. The quasar is on the 
the right hand side of the plot. The scale on the horizontal 
axis corresponds to the distance from the quasar in proper 
Mpc. The top panel shows the normalised density A along 
the line-of-sight. In the subsequent panels the output from 
the radiative transfer code is shown. In descending order 
the panels display the H i, He ii and He iii number densities 
normalised by the total hydrogen number density, the gas 
temperature, the peculiar velocity (drawn from the hydrody- 
namical simulation) , the H i, He i and He ii ionization rates 
and the energy input per ionization for H i. He i and He ii. 
The last panel shows the resulting Lya near-zone spectrum. 
The ionizing source has an age tq = 10^ yrs, luminosity 
N = 2 X 10^^ s~^ and UV background model 4 has been 
adopted. 

There are a few interesting points to note about fig- 
ure|l]before proceeding. Firstly, the H ii I-front has reached 
~ 11.8 Mpc, which is somewhat smaller than the value of 
Rion ~ 16.6 Mpc one predicts with equation ^ assuming 
A = 1. This difference is largely due to the clumpy nature 



of the IGM and to a lesser extent the presence of helium. 
Adopting A — 2.2, the actual mean normalised density along 
the line-of-sight, equation predicts Tiion = 11.8 Mpc, 
consistent with the simulation results. It is the mean den- 
sity of the gas which limits the propagation of the I-front in 
this instance. In comparison, the He iii I-front has reached 
the significantly smaller distance of ~ 4.1 Mpc from the 
source. The large difi^erence in the extent of the He iii and 
H II I-fronts is due to the smaller fraction of H i compared 
to He II for the assumed metagalactic UV background. 

The amount by which the quasar He iii I-front lags 
behind the H ii I-front depends sensitively on the H i to 
He II ratio in the surrounding IGM. Note that for a fully 
neutral IGM the extent of He iii and H ii I-fro nts should be 
very similar for typical quasar spectral index iMadau et alJ 
Il999h . There is also evidence for a clear break in the mean 
temperature of the IGM around 4 Mpc, corresponding to 
the He iii I-front position. The higher temperature is due 
to radiative transfer effects boosting the average energy 
per ionization when the He ii is reionized by the quasar 
jAbel fc H aehnelt" l999l : lBolton et al]l2004ft . This effect be- 
comes progressively larger for higher neutral helium frac- 
tions. Hence the widths of Lya absorption features ob- 
served in high reso lution spectra of quasar near-zones {e.g. 
iBecker et al.ll200^ and their correlation (or lack of) with 
distance from the quasar, combined with a lower limit to 
the extent of the H ii I-front, contain interesting informa- 
tion on the ionization state of hydrogen and helium in the 
IGM as well as the spectral index of the source. However, 
this may be complicated by the presence of high temperature 
shocked gas near the quasar host halo. Detailed modelling 
will be needed before any firm conclusions can be drawn 
from observational data. 

The second and third last panels in figure show the 
behaviour of the combined quasar and metagalactic back- 
ground radiation fields at three frequencies, corresponding 
to the ionization thresholds of H i, He i and He ii. The 
quasar radiation field falls in intensity due to increasing 
distance from the quasar and the decreasing photon mean 
free path, until the ionization rate per atom/ion drops to 
the uniform value adopted for the UV background model. 
The filtering effect of radiative transfer on the spectrum of 
the quasar radiation field is clearly visible in the second 
last panel. Higher energy photons have a longer mean free 
path, boosting the average energy per ionization in places 
where the gas is optically thick and the spectrum has con- 
sequently hardened. Note that the energy per ionization 
for H I actually decreases somewhat between the H ii and 
He III I-fronts because the photo-heating rate for H i falls 
in this region as high energy photons are soaked up by the 
optically thick helium. 

Finally, the bottom panel in figure^Jshows the resulting 
Lya spectrum for a pixel size of ~ 0.25 A. A signal-to-noise 
per pixel of 20 at the continuum level has been adopted, 
with a constant read-out signal-to-noise of 100 per pixel. 
Spectra of similar quality are already obtain able with cur- 
rent instruments jBecker et alJ 120051 l2006bl) . However, for 
the synthetic spectra used in this paper we adopt the mod- 
erate spectral resolution shown in figure m which IS more 
typical of most of the observational data. The solid vertical 
lines in both panels in figure|5|mark the extent of the quasar 
H II region at 11.8 Mpc. Clearly the extent of the Lya near- 
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Figure 4. An example line-of-sight used in our analysis of near-zone sizes. The top panel shows the normalised density distribution along 
the line-of-sight drawn from a large liydrodynamical simulation at 2: = 6 against distance from the quasar; the quasar host halo is on the 
right hand side of the diagram. Subsequent panels show the output from the radiative transfer algorithm for a quasar of age = 10^ yrs 
and ionizing luminosity N = 2 X 10^''^ s^^ . UV background model 4 was used, producing an IGM with a volume weighted neutral hydrogen 
fraction /hi ~ 10^^. In descending order the panels display the H I, He II and He III number densities normalised by the total hydrogen 
number density, the gas temperature, the peculiar velocity (drawn from the hydrodynamical simulation), the H I, He I and He II ionization 
rates and the energy input per ionization for H I, He I and He II. The last panel shows the resulting Lya spectrum with a resolution and 
signal-to-noise resembling observed high quality, high resolution quasar spectra. 
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Figure 5. The sizes of Lya near-zones around quasars at z = 6 computed from radiative transfer simulations in an inhomogeneous 
hydrogen and heUum IGM. In the left panel the IGM neutral hydrogen fraction is varied for a fixed quasar age and luminosity, in the 
central panel the quasar luminosity is varied for fixed quasar age and IGM neutral hydrogen fraction, and in the right panel the quasar 
age is varied for fixed quasar luminosity and IGM neutral hydrogen fraction. The values of the fixed quantities are listed in the bottom 
right of each panel. The sizes of Lya near-zones for individual lines-of-sight are shown with the open diamonds, with the average for each 
parameter set shown by the filled diamonds. The solid curves correspond to R™™ for = 0.4, the dashed curves correspond to /?ion 

for A = 1 and the dotted line indicates the length of the simulated sight-line. The filled grey circles and open circles in the central panel 
show the observed sizes of Lya near-zones, as determined by Fan et al. (2006b) and Wyithe et al. (2005) respectively. 



zone, Ra — 8.3 Mpc, is some way behind the H ii I-front. 
Note also that the extent of the LyP and Lya near-zones 
are the same even though the volume weighted neutral frac- 
tion of the IGM in this line-of-sight is ~ 10~^ prior to the 
quasar switching on. In this instance, naively using equa- 
tion 1121 1. one might expect the Ly/3 near-zone to be some- 
what larger than the Lya near-zone. The reason for the simi- 
larity in their size is clear on examining the density distribu- 
tion along this particular sight line, shown in the top panel 
in figureUl The transmission peaks in the Lya and Ly/3 spec- 
tra at ~ 8.3 Mpc correspond to the very low density void at 
the same position in the density distribution. Beyond this 
there are several high density clumps which attenuate the 
quasar radiation field and prevent any further transmission 
in either spectra. Consequently, variations in the density dis- 
tribution of the IGM, as well as the potential insensitivity 
of the LyQ near-zone size to the neutral hydrogen fraction, 
will make attempts to determine /hi using near-zone sizes 
difflcult. Combining the Lya and Ly/3 near-zone sizes ob- 
viously helps, but this will not necessarily give the correct 
extent of the quasar H ii region for any individual spectrum 
either. As we will discuss in more detail later, reasonably 
large samples of high quality data will be required to bet- 
ter determine the ionization state of the surrounding IGM 
using near-zone sizes. We further note that the presence of 
the lower redshift Lya forest in the Ly/3 spectrum clearly 
reduces the amplitude of the Ly/3 transition and in some 
cases will block transmission peaks from the Ly/3 near-zone. 



4.2 The measured sizes of Lya near-zones in 
cosmological density fields 

In figure|S]we explore the dependence of Lya near- zone sizes 
on the ionization state of the surrounding IGM and the lu- 
minosity and lifetime of quasars using realistic density distri- 
butions drawn from our cosmological hydrodynamical sim- 
ulation. The left panel shows the dependence of Lya near- 
zone sizes on the volume weighted IGM neutral hydrogen 
fraction. For each neutral hydrogen fraction we have per- 
formed a full radiative transfer calculation for the five dif- 
ferent lines-of-sight drawn from our simulation box at z — 6. 
The near-zone sizes were measured out to our fiducial fiux 
limit of Fiim = e~^'^ = 0.1, and are shown by the open dia- 
monds in figure Note that some of these data points fall 
on top of each other, giving the appearance of less than five 
data points for some of the parameters. The resolution of 
each synthetic spectrum was degraded and noise was added 
to mimic the observational data. The larger filled diamonds 
show the average value of Ra for each parameter set. The 
neutral hydrogen fraction of the surrounding IGM was var- 
ied using UV background models — 10, described in ta- 
ble The lifetime and luminosity of the quasar were fixed 
at tQ = lO'^ yrs and N = 2x 10^^ s"\ respectively. The solid 
curve shows 7?™'"' assuming Aii,„ — 0.4. The dashed curve 
shows iiion assuming A = 1 and the horizontal dotted line 
corresponds to the length of each simulated line-of-sight. 

The first point to note about figure|S]is the considerable 
spread in the inferred near-zone sizes from one line-of-sight 
to the next when all other parameters are fixed. This high- 



© 0000 RAS, MNRAS 000, 000-000 



12 J.S. Bolton & M.G. Haehnelt 



lights the importance of taking into account the effect of a 
reahstic density field containing voids and overdense clumps 
on the observed size and appearance of the near-zone. In 
particular, as in figure |5| it is the most underdense voids 
within the quasar H ii region which set the observed extent 
of the spectroscopically identified Lya near-zone. 

We can now investigate how the actual near-zone sizes 
relate to our analytical estimates for _Ric,n and i?™''" given 
in section |5| Firstly, we consider the effect of varying the 
surrounding IGM neutral hydrogen fraction, as shown in the 
left hand panel of figure |3 For large values of the neutral 
hydrogen fraction the observed Lya near-zone size can be 
approximated by -Rion, albeit with an offset from our chosen 
value of A = 1 due to gas clumping. The near-zone size 
increases with decreasing neutral fraction as /h/''^i but note 
also the slight deviation of the data point for an entirely 
neutral IGM from the expected scaling due to the presence of 
an extended damping wing from the Gunn-Peterson trough. 
As expected the observed sizes saturate at 7?™'"' where they 
do not increase further with decreasing neutral fraction until 
/hi < 10"^ ''. However, even for /hi > 0.1 it is difficult 
to determine the neutral fraction accurately from a small 
sample of absorption spectra due to the rather large scatter 
in near-zone sizes. Evaluating equation 11311 using i?™™ — 
7.85, corresponding to Aum = 0.4 in equation lllll . and 
/hi > 0.1 results in a timescale for the H ii I- front to reach 
Ra greater than or similar to the assumed quasar age of 
tq = 10^ yrs if A > 1. For the chosen values of the quasar 
ionizing luminosity and lifetime, saturation thus occurs at 
/hi < 0.1. For larger quasar ionizing luminosity or longer 
lifetime the saturation would occur at larger neutral fraction 
and vice versa. 

For neutral hydrogen fractions smaller than 10"'^'^, the 
Lya near-zone sizes measured using our definition (see sec- 
tion 2.2) increase rapidly as transmission from the IGM ion- 
ized by the metagalactic UV background rather than the 
quasar becomes more common. Strictly speaking, these sizes 
no longer correspond to the regions impacted by the quasar 
ionizing radiation. With further decreasing neutral fraction 
the near zone sizes rapidly approach the box-size of the nu- 
merical simulati on shown by the d otted line. This is one of 
the reasons whv lFan et alj ll2006lJ) take the necessary step 
of smoothing their Lya spectra to obtain the Lya near-zone 
sizes; not doing so makes distinguishing between the many 
thin transmission peaks from the highly ionized, underdense 
IGM and the true edge of the near-zone very difficult. How- 
ever, smoothing will have the effect of attributing transmis- 
sion peaks from the IGM to the size of the quasar near-zone, 
influencing the perceived evolution in smoothed near-zone 
sizes in the regime where the IGM makes the transition from 
being optically thick to optically thin. One must be aware of 
this efi'ect when interpreting data which has been smoothed. 
We discuss this in further detail in section |3 

In the middle and right panels of figure |S] we show the 
dependence of the Lya near-zone size on quasar ionizing 
luminosity and quasar lifetime for an IGM neutral hydrogen 
fraction /hi — 10"'' (UV background model 6). The dashed 
and solid curves again show Rion and 7?™'"', respectively. 
The simulated near-zone sizes scale with ionizing luminosity 
as Ra oc N^^'^ and are independent of the lifetime of the 
quasar as expected for this neutral hydrogen fraction, where 
Ra = -Rq Note that for a significantly neutral IGM with 



/hi > 0.1, one would see an evolution in the near-zone size 
roughly proportional to until tq > 10^ yrs. Overall, the 
analytic approximations _Rion and i?™"" appear to provide 
a reasonable guide for understanding the nature of quasar 
near-zones for 10^^ < /hi < 1. 

The filled grey circles and open circles in the middle 
panel sho w the obse r ved siz es of Lya near-zone s as d eter- 
mined by iFan et al.l J2006bl) and IWvithe et al] J2005ll re- 
spectively. There is good agreement between the smaller 
near-zone sizes and the model, suggesting consistency with 
/hi ~ 10"'' in these cases. The larger near-zone sizes, pre- 
dominantly at z < 6, lie somewhat above the values from 
the simulated spectra, but the uncertainties in the source 
luminosity, the measured near-zone size and the near-zone 
size distribution due to cosmic variance should be kept in 
mind. We will come back to this point later. Note also again 
that the simulated spectra discussed here have not been 
smoothed. As we will see in section |S] if we smooth the 
spectra this discrepancy becomes larger. 

4.3 Comparing Lya and Ly/3 near-zone sizes 

In the last section we demonstrated that Ra alone is not 
very sensitive to the ionization state of the surrounding 
IGM for typical quasar ionizing luminosities and lifetimes. 
This is due to the saturation of the near-zone size for a low 
neutral hydrogen fractions at i?™^" and the large scatter 
in near-zone sizes due to density fluctuations. The relative 
size of Ra and Rfs is an independent probe of the ioniza- 
tion state which stays useful when Ra is saturated at i?™''". 
In flgure |S] we compare the Lya and Ly/3 near- zone sizes 
predicted by 20 of the radiative transfer simulations with 
progressively smaller IGM neutral hydrogen fractions. The 
solid line shows Rp = 2.5Ra, which is the approximate up- 
per limit to the expected size of the Ly/3 near-zone pre- 
dicted by equation 11211 . The dotted line shows Ra = Rp, 
the expected relation when the IGM is substantially neutral. 
As noted earlier, these theoretical limits ignore the effect of 
foreground absorption from the lower redshift Lya forest on 
observed Ly/3 near-zone sizes. However, we do include this 
effect in the synthetic spectra analysed here. 

For /hi > 10"'^ the simulated spectra appear to fol- 
low the Rp = Ra relation. For small neutral hydrogen frac- 
tions (/hi ~ 10"'') the data exhibit an upward deviation 
approaching but not reaching the R/s — 2.5Ra relation. This 
suggests that the inhomogeneous IGM, and to a lesser ex- 
tent the additional absorption from the Lya forest in the 
Ly/3 spectrum, combined with limited spectral resolution, 
reduces the Ly/3 near-zone sizes more strongly than the 
Lya near-zone sizes. Hence, for /hi ^ 10~^ (UV background 
model 4) the Lya and Ly/3 near-zone sizes can still be rather 
similar. 

Note that for some of the completely neutral IGM mod- 
els (UV background model 0), R/s is slightly larger than 
Ra due to the Lya Gunn-Pete rson trough damping wing 
(e.g. iMesinger fc Haimanll2004l) . In the absence of this ef- 
fect one would expect Ra — Rp, since Ra has still not 
reached its maximum and thus both near-zones trace the ex- 
tent of the H II I-front. Note further that the time scale for 
the H II I-front to reach R/s is around 16 times longer than 
the timescale to reach Ra ■ Consequently, the Ly/3 near-zone 
keeps increasing in size for some time after the Lya near- 
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Figure 6. The ratio of Lya to Ly/3 near-zone sizes for in- 
homogeneous IGM models with differing volume weighted neu- 
tral hydrogen fractions. An ionizing photon production rate of 
N = 2 X 10^^ and age of tq = 10^ yrs has been assumed for 
the quasar. The solid line shows Rfj = 2.5Ra, which is the ap- 
proximate upper limit to the size of the Ly/3 near-zone predicted 
by equation I12i . while the dotted line shows Ra = Rfi- The filled 
circle with error bars show the observationally determined sizes 
of the hya and Ly/3 near-zones around J1030 + 0524 (White et 
al. 2003). 



zone has ceased to expand. This can be seen by comparing 
the results for UV background models 4 and 6. 

The filled circle with error bars show the observation- 
ally determined size of the Lya and Ly/3 near-zones around 
J1030 -I- 0524 taken from White et al. (2003). We have as- 
sumed an uncertainty in the quasar redshift of Az — 0.03. 
Its location close to the Ra — R/s line appears to suggest 
that /hi > 10~^. However, the large scatter in near-zone 
size makes it very difficult to come to definite conclusions 
for individual systems. 



5 THE OBSERVED REDSHIFT EVOLUTION 
OF NEAR-ZONE SIZES 

5.1 The simulations 

In a recent paper iFan et all i2006lJ) (hereinafter F06b) 
present evidence for evolution in the size of Lya near-zone 
sizes with redshift. The conclusion of F06b using a sam- 
ple of 16 quasars, was that the average Lya near-zone size 
increases by a factor of 2.8 from z = 6.4 to 2: = 5.7 af- 
ter correcting for quasar luminosity differences. Assuming 
that Ra oc (1 + z)^^ f^^^^ , F06b note this corresponds to 
/hi decreasing by a factor of ~ 14. Assuming a value of 
/hi ~ 9.3 X 10~^ at z = 5.7 based on their direct optical 
depth measurements suggests /hi ~ 1.3 x 10~^ at z = 6.4. 
However, as we have shown already, the dependence of the 
Lya near-zone size on /hi is not so straightforward. We shall 
attempt to model this evolution here using three different 



scenarios for the IGM neutral hydrogen fraction evolution 
within our radiative transfer implementation. 

In all three scenarios, we shall consider a quasar embed- 
ded in an inhomogeneous hydrogen and helium IGM. We set 
the quasar to have the fiducial values iV = 2 x 10^^ s~^ and 
tQ — 10^ yrs. As before, the IGM density distributions used 
for the radiative transfer runs are drawn in different direc- 
tions around the 5 most massive haloes in the simulation 
volume at z — [6.25,6,5.75], giving a total of 15 different 
lines-of-sight for comparison to the F06b data. In the first 
scenario we simulate the propagation of radiation from the 
source into the IGM at z — [6.25, 6, 5.75] using UV back- 
ground models 6, 9 and 10 respectively. These UV back- 
ground models reproduce volume weighted neutral hydrogen 
fractions which are consistent with estimates from the direct 
Gunn-Peterson optical depth measurements of F06b. This 
scenario corresponds to a moderate evolution of the IGM 
neutral fraction with redshift. In the second scenario we use 
UV background models 1, 6 and 10 at z = [6.25,6,5.75], 
respectively. This scenario corresponds to a rapid evolution 
in the IGM neutral hydrogen fraction above z = 6, where 
UV background model 1 corresponds to a substantially neu- 
tral IGM at z = 6.25, consistent with t he l ower limits on 
/hi fro m both lMesinger fc Haimaiil l|200^ and lWvithe et all 
||2005|). In the third scenario, we also adopt a moderate evo- 
lution in the neutral hydrogen fraction, but in addition take 
into account the effect of self-shielded dense clumps. These 
systems are expected to be responsible for Lyman limit sys- 
tems, but our hydrodynamical simulations do not fully re- 
solve them. To include them we use a simple model for self- 
shielded clumps base d on the model of rei onization in the 
post-overlap phase bv lMiralda-Escude et al.i ( 2000.) . The de- 
tails of this model are described in appendix D. 



5.2 The relative sizes of Lya near-zones and their 
evolution with redshift 

As discussed in previous sections and by F06b, determin- 
ing the absolute sizes of the near-zones is plagued with a 
range of uncertainties. One should, however, hope that these 
uncertainties depend weakly on redshift. Therefore, in this 
section we will firstly investigate the redshift evolution in 
the relative sizes of Lya near-zones found by F06b. We un- 
dertake a comparison of the absolute near-zone sizes with 
the observations of F06b in the next section. The evolution 
of the synthetic Lya and Ly/3 near-zone sizes with redshift 
for our three scenarios is shown in figure |7| similar to figure 
13 in F06b. We note however that figure 13 in F06b is not 
corrected for the quasar luminosity, whereas for the simula- 
tions used in figure |7| we have assumed a constant value of 
iV = 2 X 10^^ s~^. In addition, the spectra in each bin in 
figure 13 of F06b are from a range of redshifts, rather than 
all at the same redshift. 

Following F06b, the simulated spectra are smoothed to 
a resolution of 20 A and the composite absorption profile of 
the 5 lines-of-sight at each redshift is plotted. To compare 
these data to the observational results of F06b, we adopt 
their definition of the Lya near-zone. This is defined as the 
extent of the smoothed spectrum where the normalised flux 
is greater than 0.1. This hopefully provides a more robust 
way of defining Ra when we are considering examples where 
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Figure 7. The evolution of Lya and Ly/3 near-zone size with redshift for different IGM neutral hydrogen fractions. Each panel shows 
the composite absorption profile constructed from 5 lines-of-sight in redshift bins at z = [6.25,6,5.75]. The spectra have been degraded 
to resemble observational spectra and then smoothed to a resolution of 20 A. The Lya profiles are shown on the left hand side, with the 
Ly/9 profiles on the right. In all cases the quasar has an age of (q = 10^ yrs and an ionizing photon production rate of TV = 2 X 10^^ . 

Top: Near-zone sizes for an IGM ionized by a uniform metagalactic UV background, corresponding to the moderate evolution scenario. 
Centre: The UV background models are now chosen to produce a rapid evolution in the neutral hydrogen fraction. Bottom: Near-zone sizes 
for an IGM which is highly ionized at low densities by a uniform UV background but is completely neutral above a fixed density threshold. 
This mimics the behaviour of dense self-shielded clumps responsible for Lyman limit systems. 
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transmission from the thinning Lya forest affect the extent 
of the near-zones. 

The data in the top panels in figure correspond to 
the moderate evolution scenario. Between z — 5.75 and 
z = 6.25, the Lya near-zone size decreases by a factor of 
2.1. This is slightly less than the amount of evolution seen 
by F06b from z — 5.7 to 2; = 6.4, but is consistent with 
the their data over the smaller redshift range we consider. 
Using the theoretical scaling Ra oc (1 -I- z)~^ f^^^^ , this 
would imply that the IGM neutral hydrogen fraction has 
changed by a factor ~ 7 over the range 5.75 < z < 6.25. 
However, the average volume weighted IGM neutral hydro- 
gen fraction used in each redshift bin in the top panels is 
/hi = [1.6 xlO-^ 1.5x10"'', 8.1x10"^] at 2 = [6.25,6,5.75]. 
Thus, the IGM neutral fraction we use actually changes by 
a factor of ~ 20 over this range. Note that there is lit- 
tle evolution in the composite Lya near-zone sizes between 
z = 6.25 and z = 6 because the sizes of the near-zones 
are just leaving the saturated regime. A substantial part of 
the evolution of the near-zone size therefore occurred for the 
rather small change in neutral fraction between z = 5.75 and 
z = 6, where the IGM neutral hydrogen fraction changes 
by less than a factor of two. In this redshift range the 
Lya near-zones sizes are no longer near the saturated regime, 
which occurs approximately for 10~^ < /hi < 10~^. In this 
instance, transmission from the IGM ionized by the UV 
background contributes to the Lya near-zone size in the 
smoothed spectra; the flux from the quasar at the edge of 
the near-zone is comparable to the flux of the metagalactic 
UV background ionizing the surrounding IGM. Therefore, 
relatively small changes in the highly ionized IGM neutral 
hydrogen fraction can produce a rapid evolution in the ob- 
served sizes of Lya near-zones using the F06b definition of 
near-zone sizes. This occurs at a neutral hydrogen fraction of 
/hi ~ 10~*. Rapid evolution in the smoothed Lya near-zone 
sizes over the redshift range 5.75 < z < 6.25 is therefore not 
necessarily strong evidence for a significantly neutral IGM 
just above z = 6. We briefly note that we have tested this 
further by varying the flux limit used to deflne the smoothed 
near-zone sizes. For Fum > 0.1, the redshift evolution in the 
smoothed near-zone sizes is somewhat weaker. It should be 
an interesting test to perform the same exercise on the ob- 
served data. 

In the top right panel we show the corresponding com- 
posite absorption proflles for the Ly/3 near-zones. The level 
of transmission in the Ly/3 near-zone is 2 — 3 times lower 
compared to the Lya near-zone. This is mainly due to fore- 
ground Lya absorption from gas at lower redshift. The F06b 
deflnition for Ra is therefore not appropriate as a definition 
for R/s. The transmission profile of the Ly/3 near-zone is 
much shallower, more extended and fluctuates strongly. The 
Ly/3 near-zones all extend beyond the extent of the corre- 
sponding Lya near-zones, albeit by only around one Mpc 
at z — 6.25. This is expected since the optical depth for 
Ly/3 absorption is smaller and larger neutral hydrogen frac- 
tions can be probed. 

In the central panels we show the results for the rapid 
evolution scenario, with /hi = [3.1 x 10~^, 1.4 x 10~^,8.1 x 
10^^] at z — [6.25,6,5.75]. In this case, the composite 
Lya near-zone size decreases by factor of 3.1 from z — 5.75 
to 6.25. This is somewhat larger than the F06b observed evo- 
lution. Again the scaling Ra cx: (1 + z)~^ f^^^^ would imply 



that /hi has changed by a factor of 24, compared to an ac- 
tual change by a factor of over 3800. However, most of the 
near-zone size evolution in the centre-left panel has again 
occurred for the rather small change in neutral fraction be- 
tween 2 — 5.75 and z = 6. In the centre-right panel, at 2 = 6 
the Ly/3 transmission extends beyond the Lya transmission 
as before. In the 2 = 6.25 case the Ly/3 transmission extends 
the same distance as the Lya transmission, since in this case 
the edge of the near-zones correspond to the quasar H 11 I- 
front. As discussed earlier, the Lya and Ly/3 near-zones are 
expected to have a similar size for an IGM with /hi > 10~^ 
within our models. A comparison of the Lya and Ly/3 near- 
zones sizes for a large sample of high quality quasar spectra 
should thus tell whether the IGM is significantly neutral just 
above 2 = 6 or not. 

Finally, in the bottom panels we show the mean ab- 
sorption profiles for models which also take into account 
the effect of self-shielded clumps (SSCs), as described in 
appendix D. These models should be a more realistic 
representation of the post-overlap phas e of reionization 
iMiralda-Escude et alj200d : lGnedinl2000l) . The average vol- 
ume weighted neutral hydrogen fraction in each redshift bin 
for these models is /hi = [7.4x 10"^, 1.1 x 10"^, 4.3x 10"^] at 
2 — [6.25, 6, 5.75]. Note, however, that the volume weighted 
neutral fractions are considerably lower if the SSCs are ex- 
cluded, i.e. for the IGM with A < Assc. Excluding re- 
gions with A > Assc the neutral hydrogen fractions in the 
highly ionized IGM are /hi (A < Assc) = [7.9 x 10-*, 7.3 x 
10"^3.5 X 10"^]. These values are similar to those used for 
the moderate evolution scenario in the top panels. Com- 
pared to the models with no SSCs in the top panels, there 
is little difference in the Lya and Ly/3 near-zones shown in 
the bottom panels. Note that although the volume weighted 
neutral hydrogen fraction is high in these sight-lines, espe- 
cially at z — 6.25, the low density regions which dominate 
the transmission features are still highly ionized. Therefore, 
the inclusion of SSCs appears to have little effect on the 
propagation of the quasar I-front. This is perhaps not too 
surprising as the SSCs would have to contain a significant 
fraction of the total hydrogen mass in neutral form in or- 
der induce a significant evolution in the near-zone sizes with 
redshift. For a highly ionized IGM the SSCs are therefore 
not important. This may, however, be different for smaller 
neutral fractions and/or at higher redshift. 

In summary, the relative sizes of Lya near-zones and 
their evolution with redshift are consistent with a moder- 
ate change in the IGM neutral fraction. Using an analysis 
similar to that performed on observed near-zone spectra, we 
find even small changes in the neutral fraction of a highly 
ionized IGM can produce a rapid evolution in the sizes of 
Lya near-zones if the fiux from the quasar at the edge of 
the near-zone is similar to the fiux of the metagalactic UV 
background at the Lyman limit. However, because of the 
insensitivity of near-zones to the neutral fraction and the 
scatter in near-zone sizes due to fluctuations in the IGM 
density field, a substantially neutral IGM at 2 > 6 is also 
consistent with the data. A detailed analysis of the corre- 
sponding sizes of Ly/3 near-zones should help to distinguish 
between these possibilities. 
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Figure 8. The sizes of hya near-zones as a function of redsiiift. 
The filled black circles correspond to the individual Lyo near-zone 
sizes taken from the synthetic spectra in the moderate evolution 
scenario. Here the quasar is assumed to have an ionizing lumi- 
nosity of TV = 2 X 10^^ s~^. The observed correlation between 
hya near-zone size and redshift, determined empirically by Fan 
et al. (2006b), is shown with the dashed line. This fit has been 
made with observed Lyo near-zone sizes rescaled to a common 
quasar ionizing luminosity of ~ 1.9 X lO^'^ s~^, which are shown by 
the open diamonds. The simulation data lie slightly below this fit, 
although the evolutionary trend is roughly similar. The filled grey 
circles correspond to the same synthetic sight-lines, except now 
the quasar ionizing luminosity is increased to A'^ = 5 X lO^'^ s~^. 
These data are in somewhat better agreement with the fit to the 
observational data. Note that the solid vertical lines are to guide 
the eye only. In addition, only four data points are visible for 
the A'^ = 2 X 10^'' s-i data at ^ = 6.25 and z = 5.75. In these 
instances two near-zones have been measured to have the same 



5.3 The absolute sizes of Lya near-zones and 
their evolution with redshift 

We now compare the absolute sizes of the synthetic 
Lya near-zones with the observations of F06b. In figure |H1 
we plot the individual Lya near-zone sizes used in the mod- 
erate evolution scenario against redshift as filled black cir- 
cles, similar to figure 12 in F06b. The dashed line is the 
empirically determined correlation between Ra and redshift 
quoted by F06b. The dotted line is the same relation assum- 
ing the observed near-zone redshifts were overestimated by 
Az = 0.02, which is the uncertainty quoted by F06b. The 
observed Lya near-zone sizes from which this fit is derived in 
F06b, shown by the open diamonds in figure |H1 are rescaled 
to a common absolute magnitude of M1450 — —27, assuming 
A'' cx Afi45o and Ra oc iV^/^ This is an uncertain correction 
and a rescaling by N^^^ may actually be more appropri- 
ate. Using the spectral energy distribution given by equa- 
tion (O, M1450 = -27 corresponds to iV = 1.9 x 10^'' s~\ 
which is close to the fiducial value of iV = 2 x 10^*^ s~^ 
adopted for this study. We also note briefly the linear fit to 
the observational data shown in figure 12 of F06b is given by 



Ra = 8 — 9(2 — 6) Mpc, which is reproduced in our figure|H| 
This is slightly different to equation (24) in F06b, which 
they quoted erroneously (X. Fan, private communication). 

Interestingly, we obtain Lya near-zone sizes which are 
20-50 per cent smaller than the observational data. The rel- 
ative evolution and scatter in the near-zone sizes is broadly 
consistent with the F06b results. This may indicate that 
the actual density and temperature distribution of the IGM 
is somewhat different from that in our numerical simula- 
tions and that the actual mean neutral fraction is somewhat 
smaller than we have assumed for the simulated spectra. 
Another obvious way of reconciling the synthetic near-zone 
sizes of our moderate evolution scenario with the observa- 
tional data is to increase the ionizing flux in the the near- 
zone. The flUed grey circles in figure |H| correspond to the 
same synthetic sight-lines, except now the quasar ionizing 
luminosity is set to iV = 5 x lO^'^ s~^, a factor of 2.5 larger 
than before. The near-zone sizes are now in better agreement 
with the observational data. However, it appears unlikely 
that the observed quasar luminosities are underestimated 
by a factor of two to three. If an enhanced ionizing luminos- 
ity really were required to reproduce the absolute sizes of 
observed Lya near-zones, this may instead be attributed to 
star forming galaxies which are likely to cluster in the dense 
environment s urrounding the quasar. 

Recently. lYu fc Lul i200fJl suggested that the majority 
of the neutral hydrogen surrounding a quasar may have 
already been ionized by stellar sources before the quasar 
switches on. However, a factor three larger ionizing luminos- 
ity appears again too large, considering quasars are already 
extremely luminous objects. There is, however, an impor- 
tant effect which our simulations do not take into account 
which can increase the combined ionizing flux from galaxies 
(and faint quasars) in the near-zone without an enhanced 
ionizing luminosity from individual galaxies. The enhanced 
ionizing flux from the quasar in the near-zone will ionize 
dense neutral regions which otherwise block the ionizing ra- 
diation from other sources of ionizing photons. This leads 
to a significant increase of the ionizing photon mean free 
path within the near-zone compared to the mean free path 
in the surrounding IGM. This will lead to an increase in the 
metagalactic UV background flux due to the galaxies and 
possibly other, fainter quasars within the near-zone. An ac- 
curate estimate of this effect is difficul t but a rough estimate 
can be obtained using our SS C model iMiralda-Escude et alJ 
I2OOOI: iFurlanetto fc Ohll2005l and also see appendix D). An 
increase by a factor of two or three in the ionizing fiux at 
the edge of the near-zone appears plausible. 

There is a range of other effects which may affect the 
estimated near-zone sizes. However, they should be less im- 
portant than the effects discussed above. Spatial fluctua- 
tions in the UV background expected around the end of 
reionization may influence the size of observed near-zones. 
In addition, gas in highly overdense regions is likely to be 
shocked, producing hot, highly ionized gas. This has not 
been taken account in our simulations. However, the vol- 
ume fllling factor of this hot, ionized gas phase is expected 
to be small at z — 6. The difficulty of the data reduction 
may also contribute to the discrepancy. The continuum fit- 
ting of these highly absorbed spectra in the region of the 
broad Lya emission line is notoriously difficult. Unidenti- 
fied emission lines could easily lead to a quasar continuum 
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which is placed too low at the relevant distance from the sys- 
temic redshift. This would result in an overestimate of the 
transmission and could therefore result in an overestimate 
of the observed near-zone sizes. 



6 CONCLUSIONS AND DISCUSSION 

We have presented a detailed investigation of the properties 
of near-zones in the Lyman series absorption spectra of 
z ~ 6 quasars, using analytical estimates and an accurate 
radiative transfer scheme for a hydrogen and helium IGM 
combined with density distributions drawn from a largo 
hydrodynamical simulation of a ACDM universe. From the 
radiative transfer simulations we have produced realistic 
mock spectra which we have analysed in a similar manner 
to the current observational data. The main conclusions of 
this study are as follows. 

• For large IGM neutral hydrogen fractions the size of 
the near-zone increases eis /h/^^, as expected for an 
H II I-front expanding into a significantly neutral medium. 
For typical luminosities of observed z ~ 6 quasars and 
plausible lifetimes, the near-zone sizes saturate at a neutral 
fraction /hi w 0.1. For smaller neutral fractions the spec- 
troscopically identified near-zones are more akin to classical 
proximity zones within a highly ionized IGM, and their size 
is largely independent of the neutral hydrogen fraction of 
the surrounding IGM. 

• The density fluctuations in a cosmologically representative 
matter distribution lead to considerable scatter in the size 
of the near-zones, even for fixed assumptions about the 
ionization state of the surrounding IGM. It is thus difficult 
to determine the IGM neutral hydrogen fraction from a 
small number of observed spectra. 

• Observed sizes of spectroscopically identified near-zones 
in the Lya region alone axe a rather poor indicator of 
the ionization state of the surrounding IGM. Observed 

near-zone sizes can furthermore substantially underestimate 
the size of the region impacted by the ionizing radiation of 
the quasar. 

• For small neutral hydrogen fractions quasar He iii I-fronts 
lag considerably behind their H ii counterparts if He ii is 
still abundant. This leads to a temperature jump behind 
the He iii I-front which can be located within the observed 
near-zone for a suitable combination of IGM neutral frac- 
tion, quasar luminosity and quasar lifetime. Detailed study 
of the width of the absorption features in the Lya near-zone 
may therefore reveal the position of the He iii front. This 
could then be used together with the Lya near-zone size to 
constrain the neutral hydrogen fraction. 

• The ratio of near-zone sizes in the Lya and Ly/3 region 
of quasar spectra provides another independent probe of 
the ionization state of the surrounding IGM. For neutral 
fractions /hi > 10^^ we find that the size of the Lya and 
Ly/3 near-zones are similar. For smaller neutral fractions 
the size of the Ly/3 region becomes substantially larger, 
approaching the theoretical limit of Rp ~ 2.5Ra- The 



extent of the Ly/3 near-zone relative to Lya is blurred by 
the overlying Lya forest at lower redshift. It is also very 
sensitive to the density distribution along the line-of-sight. 
For this reason the ratio of the sizes of the Lya and 
Ly/3 near-zones approaches the saturated value at a neu- 
tral fraction somewhat smaller (/hi^10~^) than where 
saturation sets in for the sizes of Lya near-zones alone 
(/hi ^ 10~^). The exact value of the threshold neutral 
hydrogen fractions depends on the quasar luminosity and 
lifetime. 

• The synthetic spectra produced from our radiative 

transfer simulations exhibit Lya near-zones which are 
20 — 50 per cent smaller than the observed spectra for 
plausible assumptions for the ionizing luminosity of the the 
quasar and the flux of the metagalactic UV background in 
the surrounding IGM. The discrepancy is resolved if the 
ionizing flux at the edge of the near-zone is enhanced by a 
factor of two to three. Such an enhancement is plausible 
given the expected increase in the ionizing photon mean free 
path in the near-zone compared to the mean free path in 
the surrounding IGM. This would produce a corresponding 
enhancement of the metagalactic UV background within 
the near-zone from stellar emission and possibly fainter 
quasars. Alternatively, it may suggest that the IGM has 
a lower neutral fraction than we assumed in our simulations. 

• The strong evolution of the observed sizes of Lya near- 
zones over the short redshift interval 5.7 < z < 6.4 is 
consistent with, but does not require, a rapid evolution of 
the neutral hydrogen fraction of the surrounding IGM. The 
observed evolution in the Lya near-zone size by around 
a factor of two is equally consistent with a moderate 
evolution in the IGM neutral hydrogen fraction. Relatively 
small changes in the highly ionized IGM neutral hydrogen 
fraction can produce a rapid evolution in the observed sizes 
of Lya near-zones if the ionizing flux from the quasar at 
the edge of the near-zone is comparable to the flux of the 
metagalactic UV background at the Lyman limit in the 
surrounding IGM. In this case the rapid evolution can be 
attributed to additional transmission from the thinning 
Lya forest. This proposition can be tested by varying the 
flux level used to define the extent of the near-zone sizes. 
For larger flux levels the evolution should decrease if the 
increasing transmission from the general Lya forest is 
responsible for part or most of the evolution. 

• The effect of self-shielded clumps on the observed sizes 

of the Lya near-zone is small, and the expected evolution 
of the incidence rate of self-shielded clumps appears to 
contribute little to the rapid evolution of near-zone sizes. 

We note here that we have not considered two further pos- 
sibilities which may influence the size of quasar near-zones: 
a fluctuating UV background and additional sources of ion- 
izing radiation, such as star forming galaxies, clustered in 
the vicinity of the quasar host halo. Both of these are likely 
to increase the observed near-zone sizes. Additionally, ob- 
servational constraints on near-zone sizes are still somewhat 
uncertain, as are quasar ionizing luminosities and lifetimes. 
However, a large sample of high quality quasar spectra 
with an accurate determination of the sizes of the Lya and 
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Ly/3 near-zones should improve the constraints on the ion- 
ization state of the IGM at z > 6 considerably. If most 
quasars exhibit Ly/3 near-zones which are similar in size to 
the counterpart Lya near-zone, this should favour an IGM 
with /hi > 10"^ If Ly/3 near-zones were instead generally 
larger than their Lya counterparts this would provide evi- 
dence for an IGM with /hi < 10~^. The information locked 
within the widths of absorption features in the Lya near- 
zones may provide an interesting additional probe of the 
ionization state of the IGM at high redshift. 
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APPENDIX A: TIME DELAY EFFECTS ON 
OBSERVED SIZES OF COSMOLOGICAL 
H II REGIONS 

The existence of a finite speed for the transmission of real 
information alters the observed sizes and shapes of cosmolog- 
ical H II regions. These time delay effects, although having 
their basis in special relativity with the postulate that the 
speed of light in empty space always has the same value, are 
purely geometrical rather than relativistic phenomena which 
manifest themselves when the I-fronts travel at relativistic 
speeds. The renewed interest in cosmological H ii regions as 
a probe of the epoch of reionization has meant these effects 
have been well documented within the r ecent literature, to 
which we refer the intereste d reader (e.g. [ White et alj|2003l : 
IYu fc Lui 2005: Yu 2005: Shapiro et al J l2006ll However, we 
briefly review these results here, as they provide an impor- 
tant simplification to our radiative transfer scheme. 

One can easily show the expression for the radius, R, of 
a H II region expanding into a static, pure h ydrogen medium 
of constant density and te mperature is (e.Q. IWvithe fc Loebl 
l2004l:IShapiro et al-llioo^ 



where A*' is the number of ionizing photons emitted by the 
monochromatic source per unit time, nn is the number den- 
sity of hydrogen atoms, /hi = JiHi/nn is the neutral hy- 
drogen fraction and qhii is the recombination coefficient for 
ionized hydrogen. This expression explicitly takes into ac- 
count the time delay for a photon to travel from the source 
to the edge of the H ii region, thus setting the limiting ex- 
pansion speed of the H ii region to the speed of light. 

Now consider the expansion of a H ii region outwards 
from the ionizing source along the line of sight to an ob- 
server. Any photons instantaneously detected by the ob- 
server were emitted at the same retarded time, = t — R/c, 
where t is the time when the photons reach a distance R 
from the source. Evaluating equation lAll at the retarded 
time yields 



dR _ N - ^■KR^aminj 
dfa ATvR'^fmnH 



(A2) 



dR 
'dt 



N -f 4nR'^ fmriHC — ■^TvR^amin'^ 



(Al) 



which is identical to equation Q. Thus, it is clear for 
this particular geometrical configuration that the expansion 
speed of the H ii region in the observer's frame is no longer 
limited to the speed of light; it is as if one can completely ig- 
nore time delay effects and assume real information is trans- 
mitted instantaneously. This leads to the observation of an 
initial superluminal expansion phase close to the ionizing 
source, when the number of ionizing photons per hydrogen 
atom at the Lfront greatly exceeds unity, allowing it to prop- 
agate at relativistic speeds. This result is particularly useful 
in the context of our radiative transfer simulations, which 
we discuss in the next section. 



APPENDIX B: RADIATIVE TRANSFER 
IMPLEMENTATION 

Modelling the extent of cosmological H ii regions around 
high redshift quasars, while consistently tracking local 
changes in the ionization state and temperature of the gas 
due to the filtering of ionizing radiation, requires a full treat- 
ment of radiative transfer through a hydrogen and helium 
IGM. In a cosmologically expanding medium I-fronts from 
a luminous quasar are highly supersonic and thus decoupled 
from the hydrodynamical response of the gas. This allows 
one to ignore the radiation hydrodynamics of this particular 
problem, providing a significant simplification to our nu- 
merical scheme. For this purpose we use an updated version 
of the radiative transfer algorithm of Bolton ct al. (200^. 
This algorithm is based on the explicit photon conserva- 
tion method of lAbel et al.l 1^9^, which has the desirable 
property that photons are conserved almost independently 
of spatial resolution, allowing the accurate determination of 
physically relevant quantities with only moderate compu- 
tational expense. We briefiy outline the updated algorithm 
here and then present some tests of the code. 



Bl The cosmological model 

We adopt the following time-redshift relation for a flat 
Universe dominate d by a cosmo logical constant and non- 
relativistic matter JPeebleslll99'3) 
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Hot 



3(1 -QmV/^ 



sinh 



V fin, / 



1/2 



3/2 



(Bl) 



where Ho ~ WOh km Mpc~^ is the present day Hubble 
parameter, ilm is the matter density as a fraction of the 
critical density and a = (l + 2)~^is the cosmological scale 
factor. We use this relation to update redshift dependent 
quantities at each successive timestep 5t within the radiative 
transfer implementation. Additionally, from equation IBlll 
it is straightforward to derive an expression for the Hubble 
parameter 



H{t) = n^/^ Ho cosh 



3Ho 



(1 - n,. 



-3/2 



(B2) 



B2 The ionizing radiation field 

Consider a source which emits ionizing photons per unit 
time, such that 



N 



Nydv ■ 



where dv is the luminosity of the source within the fre- 
quency interval v ^ u + du, defined as 



Lv du — Li 



dv. 



The source has a power law spectrum with index as, typical 
of quasars, normalised by Lhi at the H i ionization threshold 
frequency, //hi- 

The probabilities of an ionizing photon emitted by the 
source the n being absorbed b y H i, He i and He ii, respec- 
tively, are jBolton et al.ll2004^ 

-Phi = Phi(7hcI(7hcIi(1 — e " )/D, 

-PhoI = gHlPHcigHoIl(l — e " )/D, 

-PhoII = (jHigHclPHcIl(l - e " )/D, 



where pi = 1 — e , iji = e , r^°* = r^" + t;/" + 

is the optical depth for species i and D = pHi^HoigHcii + 

9HlPHeigHcII + ?Hl9HclPHeII. 

If we now invoke a uniformly discretised space grid with 
elements of size SR along the line-of-sight to the source, at 
any given timestep St, the photo-ionization rate per unit 
volume [s~^ cm~'^] of species i in grid element I is 



l-nl 



Nl'^Pi dv, 



where 1/' is the volume of element I and Ui is the ionization 
threshold frequency for species i. If we assume the luminosity 
of the source is constant and allow an additional contribu- 
tion to the ionization rate from a diffuse UV background, 
the above equation can be rewritten as 



l-nl 



1 



k 



ly^b 

- nil i , 



(B3) 



which is the equation we shall solve for the ionization rate 
per unit volume for species i. Here we have optimally dis- 
cretised the frequency integration into 20 logarithmically 
spaced intervals k between v\ and lO/^i, R is the distance 



of element I from the source, V\ is the optically thin photo- 
ionization rate due to the diffuse metagalactic UV back- 
ground and the transmission factor, Tu, is defined as 



= 0), 

I]!=o('^Hin^j -f (THoin^^i + 0-HoIin^3„)(5i? {I > 0). 



Similarly, the photo-heating rate per unit volume 
[erg s~^ cm^^] for species i is given by 



1 1 



k 



Pi/lp(fk 



+ nV\, (B4) 



such that the total photo-heating rate per unit volume in 
element I is given by = riHiem + '^HcieHoi + nHoiieHdi 
and is the optically thin photo-heating rate due to the 
diffuse UV background. 

In general, for each timestep &t the number densities ni 
in the expression for the transmission factor are evaluated 
in each grid cell j at the same retarded time tR = i — j5R/c, 
where t is the time it takes a light signal to reach cell j. In 
principle, an I-front propagating outwards from the source is 
then limited to travel no faster than the speed of light since 
rii are only evaluated when t > jSR/c. However, in prac- 
tice spatial discretisation means this is not always the case; 
an I-front is assumed to propagate instantaneously across 
a cell width SR. Even for relatively high spatial resolution, 
this enables the expansion rate of the ionized region to ex- 
ceed c. This means that to accurately compute the position 
of an I-front one must either increase the spatial resolution 
to a prohibitively high level or artificially impose the limit 
R < ct. At first glance, the latter option would seem to be 
preferable. However, this has the undesirable effect of vio- 
lating photon conservation, so that alth ough one correctly 
predicts the position of the I-front (e.g. lAbel et al. 199iJ, 
one must again resort to substantially higher spatial resolu- 
tion (such that each cell in the space grid is optically thin) 
to compute physically relevant quantities within the region 
where R < ct is imposed. This is a problem for simulat- 
ing quasar near-zones in Lyman series absorption, since the 
absorption spectrum is sensitive to the exact values of the 
neutral fraction and temperature of the gas near the source. 
Fortunately, we have already shown that in the frame of 
an observer the expansion rate of the I-front follows equa- 
tion 1A2II . Since we are interested in the sizes of observed 
H II regions, we may evaluate the expression for assuming 
the photons are transmitted instantaneously, removing the 
need to impose the condition R < ct while at the same time 
explicitly conserving photons. This enables us to determine 
the ionized fraction and temperature of the gas close to the 
ionizing source accurately for relatively coarse spatial reso- 
lution, substantially reducing the computational expense of 
the simulations. 



B3 Computing the ionized fraction and 
temperature 

At every timestep St, in each spatial grid element SR, the 
ionization state and temperature of atomic hydrogen and 
helium must be computed by solving a set of four coupled 
first order ordinary differential equations. The first three 
determine the abundances of the three ionized species 
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dt 
driHcii 



dt 



driY 



— nHl(rHI + ncrcHi) — nHIinoQ!HII, 

= «-Hci(rHci + neTcHci) + ?^HcI^?^cQHoI^ 

— T^iHcIl(rHcII + ".crcHcII + ?^oQHcIi), 

- = MHcIl(rHcII + "-crcHcIl) — riUcluncCMHc] 



dt 

with the closing conditions 



riHi — nu — nmi, 



(B5) 

(B6) 
(B7) 



Y 



riHci = 



-riH — riHcii — nHciii, 



4(1 ~y) 

ria = nmi + nuaii + 2nHGin- 

Here ni denotes the number density [cm~^] of species i, Fi 
[s~^] is the photo-ionization rate, Fci [cm^ s~^] is the col- 
lisional ionization rate, Ofi [cm^ s~^] is the total radiative 
recombination coefficient and Y is the helium mass fraction , 
which we take to be Y = 0.24 (e.Q. IOlive fc Skillmanll2004^ . 

Additionally, the temperature must also be computed 
by solving a fourth differential equation 



dT 
HI 



(7 - l)/xmH 



A{ni,T)]-2Hit)T, 



(B8) 



where fi is the mean molecular weight of the gas, 7 — 5/3, 
A(ni,r) [erg s~^ cm~^] is the radiative cooling function, 
H{t) is the Hubble parameter given by equation IIB2t and 
all other symbols have their usual meanings. The last term 
on the right hand side of equation IB8II is the contribution 
to cooling by the adiabatic expansion of the Universe. 

We solve equations <B5ll - IIB8ll using a first order implicit 
integration scheme which provides a good balance between 
accuracy and speed iAnninos et aP llQQ'iD . We find that to 
solve these equations accurately, one must typically adopt a 
timestep comparable to the hydrogen ionization timescale, 
tion ~ I/Fhi, which can become very small in close proxim- 
ity to a luminous source. On occasion this time interval is 
smaller than the main timestep we use to track the position 
of the I- front (see section [HTI for details). When this is the 
case we sub-cycle over each main time step using smaller 
time intervals based on the ionization timescale to solve for 
the ionic abundances. Once the rate of change in the elec- 
tron number density is less than 10~^, we switch to solving 
the ionic abundances assuming the gas is in ionization equi- 
librium {e.g. Katz ct al. 1996). This leaves only the differen- 
tial equation for the temperature to be integrated, allowing 
much larger timesteps to be taken which greatly speeds up 
the calculation. 



behind the source I-fronts and is thus optically thin. In ad- 
dition, for typical gas densities the recombination timescale 
is considerably longer than the canonical quasar lifetime. 
In these circumstances we find case A recombination coeffi- 
cients are the appropriate choice. 

(1) Photo-ionization cross-sections [cm ] (IOsterbrockll989ll : 



6.30 X 10"^^[1.34(i//z/Hi)" 



-18 [ 

-IBr 



0.34(i//^^i) 



o-Hci = 7.03 X 10-"°[1.66(i//i^oi)" - 0.66(i^///Hci) 



-2.05 
-2.99 



0-Hcii = 1.50 X 10-"''[1.34(i//i^o„)-^'^^ - 0.34(i//z/Hcii 
(2) Optically thin photo-ionization rates [s~'^]: 
rl^i = 1.27 X 10-"j^i(ab-f3)-\ 
rl^ei = 1.51 X 10"" (0.553)°" J^! 

T^b o no 1 n-12 



\-3.99i 



3.03 X 10" 



(ab + 2)-\ 
(0.250)""J^i(ab-f 3)-\ 



where the specific intensity of the metagalactic radiation 
field [erg s"'^ cm~^sr~^Hz~'^] is defined as 

J^ = J^iXlO-^^(^)"'^\ 

The radiation field is assumed to have a power law spec- 
trum with index ah, normalised by J^2i at the HI ionization 
threshold. 

(3) Optically thin photo-heating rates [erg s~^]: 

e^i = 2.94 X 10-''j^2i[(ab + 2)(Qb + 

eLi = 6.48 X 10-''(0.553)"V!!2i[(Qb + l)(ab + 2)]-\ 

eLii = 2.80 X 10-''(0.250)"V^2i[(ab + 2)(ab + 3)]-\ 

(2) Recombination rates [cm s' (lAbeTeTiDllii^): 

ami = exp(-28.6130338 - 0.72411256 In(f') 

2.38086188 x 10"^ln(T)^ 
1.42150291 X 10""ln(f)'^ 
-H4.98910892 x 10"'^ In(f')® 5.75561414 x 10"''ln(f)'' 



-2.02604473 x 10"^ln(T)^ 
-3.21260521 X 10"''ln(f)'' 



-1.85676704 x 10"*ln(f)^ - 3.07113524 x ln{ff , 

_ T . d 
— <^HcII + CtHcIIi 



3.925 X io-i3j.-o.6353 



r 

d 1 r-AA „ ^-90^-1.5 -48.596/f /n o , 8.1/T 

QhgII = 1.544 X 10 J e ' (0.3 -\- e ' 

QHoiii = 2aHii(f'/4), 

where T is the gas temperature in electron volts. 



(4) CoUisional ionization rates \cw? s '^l 



ijTheuns et alJ 



B4 Cross-sections and rates 

The cross-sections and rates adopted in our radiative trans- 
fer scheme are listed below. Note that we use case A instead 
of case B recombination coefficients for this study. For op- 
tically thick hydrogen and helium, one can approximately 
take into account the moderate boost to the photo-ionization 
rates due to diffuse radiation emitted by recombining ions 
using case B recombination coefficients. T his is the so-called 
on-the-spot approximation lIsiitieJigTi). However, any gas 
seen in Lya absorption in our simulated spectra typically lies 



ri 1 TT^, in-10T^0.5 -157809. 1/T r, . rrnl/2T-l 

ToHi = 1.17 X 10 T e ' [I-I-T5 J , 

r A 'vr An-llrrO.5 -285335. 4/T r, , rTil/2i-l 

roHci = 4.76 X 10 T e ' [l-fTg' J , 



roHdi = 1.14 X io-''r"-e 

where T5 = T/10^ K 



Ilrp0.5„-631515.0/T 



[1 + T^JT\ 



The following radiative cooling processes contribute towards 
the cooling function A(ni,T): 

(5) Recombination cooling [erg cm"'']: 
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Table CI. Radiative transfer test models. For each model the 
columns list the ionizing photon production rate of the source A*^, 
the spectral index of the source Cfg , the hydrogen number density 
nji, the mass fraction of helium Y and the age of the source tg. 
Models C and D are for a monochromatic source, therefore no 
source spectral indices are listed. 



Model 


N [s-i] 


as 


nil [cm ^] 


Y 


is [yrs] 


A 


5 X 10^^ 


1.5 


7 X 10-5 


0.24 


lO'' 


B 


5 X 10^"^ 


1.5 


7 X 10-2 


0.24 


10^ 


C 


1 X 10^5 




7 X 10-5 





107 


D 


1 X 10^9 




7 X 10-2 





10^ 



A'hTi = 1.036 X lO^^'^raHimenHii, 

Ahoii = (1-036 X lO^'^'TaJirfi + 6.526 x 10""aHcii)nenHoii, 
AhoIii = 1-036 x 10~^''TQHciimonHcin- 

1 6) Collisional i onization cooling [erg s-"*^ cm-'*] 
Theuns et"aIlll998^ : 

AoHI = 2.18 X 10""reHmenHI, 
AcHcI = 3.94 X lO^^^rcHcmcnHcI, 
AcHcII = 8.72 X lO^^^rcHcimcnHoII- 

(7) Collisional excitation cooling [erg s"^ cm"^] iCeDll992l) : 



A.^1 = 7.50 X lO-^-'e 



19 -118348/T 



[l + Tg/^l'^nenm 



AS^i = 9.10 X io-27r'°-^««^e-"^^«-°/^[i + r^/^i-^-^ 



71(3 nHoI, 



AhoII = 5.54 X 10 1 8 ' [l + ^^s I "-cnHoii- 



(8) Bremsstrahlung [erg s"^ cm"^] llCeJl992h : 

Afi = 1.43 X 10"^'^T^''^gfi(nHn + nncu + 4nHoni)no, 

where the Gaunt factor for free-free emission, gg, is given 
by 

ga = 1.1 + 0.34 exp[-(5.5-logio(r)]^ 

(9) Inverse Compton cooling [erg s"^ cm" JPeebleslll973) : 
Ac = 5.65 X 10"^''[r - 2.73(1 + z)]{l + zfn^. 

APPENDIX C: CODE TESTS 

Exploring the limitations of any numerical model is essential 
for providing well motivated results. To this end we present 
some selected tests of the radiative transfer implementation 
we use for our simulations. The test models are summarised 
in tableau 

CI Spatial and temporal resolution requirements 

We firstly present some convergence tests for the spatial 
and temporal resolution of our radiative transfer scheme. 
Our algorithm has the useful property that photons are 
conserved explic i tly independent of s patial resolution (e. 
lAbel et alJ Il999l: iBolton et all 12004 iMeUema et"ai]|200i 



Consequently, the ionized fraction and temperature are com- 
puted accurately even when the cells on the space grid are 



very optically thick, speeding the calculation up consider- 
ably. 

The optical depth to neutral hydrogen in each cell of 
width &R on the space grid is given by irni = cthiMhi^^?- 
Employing fiducial values, this can be written as 



5rHi = 13.6 



.7 X 10- 



10 kpc 



(CI) 



Note that in the case of a hydrogen and helium IGM, the 
optical depth to He i and He ii per cell is less than that for 
H I due to the lower number density of these species and, in 
the case of He ii, a smaller photo-ionization cross section. In 
this case the H i optical depth per cell is the correct proxy 
to use for the spatial resolution. 

The results of a spatial resolution test using models A 
and B are shown in figure IHTI We have adopted the same 
time step size,e = 0.05, in both of these models. We shall 
discuss this in more detail shortly. In both models, the prop- 
agation of ionizing radiation into a hydrogen and helium 
medium of uniform density is computed for a source with a 
power-law spectrum situated 2 kpc from the edge of the sim- 
ulation grid. The left panels show the fractional abundances 
of H I, He II and the gas temperature for model A, with a 
hydrogen number density of 7 x 10-^ cm"^, comparable to 
the mean density of hydrogen ai z — 6. Even adopting cell 
sizes in which the hydrogen is initially very optically thick 
we find all quantities are recovered with a good degree of 
convergence. The right hand panels show the same results 
for model B, which has a hydrogen number density 1000 
times greater than model A. Again, the degree of conver- 
gence for such optically thick cells is quite remarkable. Note 
there is evidence that the gas temperature has not fully con- 
verged at the lowest spatial resolutions, especially close to 
the source and near the He iii ionization front where there 
are rapid changes in the energy input per ionization. How- 
ever, we find a spatial resolution of SR ~ 10 kpc provides 
a good degree of accuracy for the range of parameters we 
shall consider. 

The second factor we need to consider is the temporal 
resolution. As already discussed in the previous section, po- 
tentially the shortest timescale we need to consider is the 
hydrogen ionization timescale tion ~ I/Fhi. However, this 
timescale varies locally in the simulation. More generally the 
timestep size we are interested in is the cell crossing time for 
real information, which we parameterise using the Courant 
condition e, defined as 



cSt 

Jr' 



(C2) 



We adopt this as the main timestep size of our radiative 
transfer scheme. Note here that the time resolution is de- 
pendent on the spatial resolution within our implementa- 
tion. Again, using typical values this can be rewritten as 



St = 3261.6 — 
.0.1 



SR 
10 kpc 



yrs. 



(C3) 



Figure IC2I shows the results of our temporal resolution 
tests for models A and B. We have adopted a cell size of 
SR — 1 kpc for all of these test runs. It is immediately 
clear from model A that a timestep equal to or shorter than 
the cell crossing time (e < 1) is required to achieve nu- 
merical convergence. For model B which has a higher gas 
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Figure CI. Test of spatial convergence for our radiative transfer scheme. The left hand panels show the fractional abundances of H I, 
He II and the gas temperature for model A using the cell sizes labelled on the plot. The right hand panels show the same quantities for 
model B. 



density, a stricter requirement e < 0.1 is required to pro- 
duce reasonably converged results. This is because the rate 
of change in the optical depth of each cell is greater per 
timestep. For fully three-dimensional cosmological radiative 
transfer with multiple sources this timestep size would be 
somewhat prohibitive and could be improved upon by using 
a for m of adaptive timest epping based on the I-front ve- 
locitv ij Shjjjiro et_al^ or time averaged optical depths 

l)Mellem^^^dTl200m . However, for our purpose this simple 
time stepping scheme is perfectly adequate. Consequently, 
we adopt a value of e = 0.1 for our simulations. 

C2 Pure hydrogen isothermal H ii region 

Modelling the expansion of a H ii region into neut ral hy- 
drogen gas of uniform density (e.g. IStromgren|[l939tl is one 

of the few radiative transfer problems which has an exact 



analytic al solution with w hich one can compare simulation 
output illiev et al.ll2006'ah . We consider two test cases (C 
and D) which simulate the propagation of an H ii I-front 
into an initially neutral pure hydrogen medium of uniform 
density and temperature (T = 1.5 x 10'' K). The ionizing 
source is assumed to be steady and monochromatic, emit- 
ting A'^ photons per unit time at the Lyman limit. 

We compare the resuhs of test models C and D to the 
classical analytical solution, equation (|3J, in figure IHsl The 
solid lines show the numerical solution for the H i fraction 
at t/lO'^ yrs = [0.2,0.4,0.6,0.8,1.0] for model C and 
t/lO'' yrs = [0.05,0.1,0.2,0.3,1.0] for model D. The solid 
grey strips give the analytical solution with a width equal 
to the cell size of the simulation, and hence the uncertainty 
in the numerical result. For this test we adopt a cell size 
of SR = 20 kpc, with e = 0.1 for model C and e = 0.01 for 
model D. The agreement between the numerical and analyt- 
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Figure C2. Test of temporal convergence for our radiative transfer scheme. The left hand panels show the fractional abundances of 
H I, He II and the gas temperature for model A using the Courant conditions labelled on the plot. The right hand panels show the same 
quantities for model B. 



ical results is very good, which suggests that we are calculat- 
ing the extent of the H ii regions correctly within our radia- 
tive transfer implementation. Note that by t/W^ yrs =1.0 
model D has reached its Stromgren radius. 



APPENDIX D: A SIMPLE MODEL FOR 
SELF-SHIELDED REGIONS 

After the low density gas in the IGM has been reion- 
ized, dense, self-shielded clumps (SSCs) are thought 
to s et the mean free path for ionizing photons 
(e.g. iMiralda-Escude et alJ l2000l . hereinafter MHROO, 
iMiralda-Escudj l2003l) ." We adopt a simple model to in- 
clude these SSCs with in our simulations. As noted by 
iFurlanetto fc OhI feOOST) . one may approximate the density 
at which hydrogen in a dense clump becomes self-shielded. 



For a clump of normalised density A 
equilibrium wit h the UV bac kground, 
approximately llSchavell200ih . 



~ pb/Pb in ionization 
its column density is 



where F^, 



s ^ is the normalised metagalactic 



hydrogen ionization rate and T is the gas temperature. The 
clump will become optically thick to ionizing radiation once 
Nm > l/o"Hi ~ 1.6 X 10^^ cm~^. Adopting this value for 
A^'hi and rearranging yields the expected normalised density 
threshold for a SSC immersed in the metagalactic radiation 
field 



49.5 



2/3 



(D2) 
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Figure C3. Comparison of test models C (upper panel) and D 
(lower panel) to the analytical solution for the size of a H II region 
in a pure hydrogen medium of constant temperature and density. 
The solid lines show the numerical solution for the H I fraction 
at t/lO"^ yrs = [0.2, 0.4, 0.6, 0.8, 1.0] for model C and t/lO"^ yrs = 
[0.05,0.1,0.2,0.3,1.0] for model D. The centre of the solid grey 
strips give the analytical solution with the width of the strips 
representing the cell size of the simulation. 



In comparison, the number of clumps per unit redshift with 
A > [4.1,21.2,37.5] at z = [6.25,6,5.75] within our 400^ 
simulation volume is around a factor of two less than this. 
Therefore, when setting all clumps with A > Assc to be 
neutral within our simulations using UVB models 6, 9 and 
10, we underpredict the mean free path for ionizing pho- 
tons by about a factor of two. To compensate for this, we 
lower the density threshold at which a clump becomes self- 
shielded, thus reproducing the required number of SSCs per 
unit redshift given by the MHROO model, as listed in equa- 
tion llD4ll . The relevant density thresholds in our 400'^ simu- 
lation which reproduce the required number of SSCs per unit 
redshift are Assc = [2.7,10.0,20.0] at z = [6.25,6,5.75]. 
Thus, our simulations using Assc = [2.7, 10.0, 20.0] and UV 
background models 6, 9 and 10 at z = [6.25,6,5.75] should 
be equivalent to using Assc ~ [4.1,21.2,37.5] with Amfp = 
[1.09,6.99,12.80] Mpc and T^ia = [0.028,0.280,0.560]. As- 
suming ionization equilibrium at A = 1 and T — 2 x 10* K, 
these ionization rates correspond to /hi = [7.3 x 10~*, 7.3 x 
10"^ 3.6x10"^]. However, we note that simulations incorpo- 
rating radiative transfer which correctly resolve the number 
density of SSCs are needed to fully address their effect on 
observed near-zone sizes. 



For UV background models 6, 9 and 10, which produce IGM 
neutral fractions consistent with the observational data of 
iFan et al.l ll2006bi) . this implies Assc = [4.1,21.2,37.5] at 
z — [6.25,6,5.75]. Ideally, we would implement this self- 
shielding for A > Assc within our radiative transfer simu- 
lations in the same way as MHROO, by assuming all regions 
with A > Assc remain neutral while still optically thick to 
ionizing radiation. However, although we have shown that 
the resolution of our hydrodynamical simulation should be 
sufficient for correctly modelling the propagation of ioniz- 
ing radiation, it still does not quite reach the required mass 
resolution to reproduce the correct number of clumps with 
A > Assc- MHROO overcome this problem by using an IGM 
density distribution drawn from a hydrodynamical simula- 
tion, with the addition of a power-law high density tail to 
account for the lack of high density clumps. The mean free 
path for ionizing photons in their model is given by 

A^fp = [1 - Fv(A < Assc)]"'^' Mpc, (D3) 

where H{z) is the Hubble parameter and -FV(A < Assc) is 
the fraction of the IGM volume with A < Assc, computed 
from the volume weighted density distribution corrected for 
high densities. Assuming Amfp is similar to the mean distance 
between SSCs, which is exp ected to be a reasonable approx- 
imation (iFurlanetto fc Obl | 2"0 05'l. the number of SSCs per 
unit redshift with A > [4.1,21.2,37.5] at z = [6.25,6,5.75] 
using the MHROO model is 

- T^TTTTT-. T = [52.5, 9.0, 5.1]. (D4) 

dz X^ip H{z){l + z) L , > J \ J 
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